#packages we need for this code file
library(ggplot2)
library(mgcv)
library(lubridate)
library(zoo)
library(tidyverse)
library(dplyr)
library(DHARMa)
library(mgcViz)
library(extrafont)
library(arm)
loadfonts()
library(stargazer)
library(ellipse)
library(dotwhisker)
library(countreg)
#define functions we will need for analysis
#expit function
expit<-function(x){
return(exp(x)/(1 + exp(x)))
}
#logit function
logit<-function(x){
return(log(x/(1 - x)))
}
#read in data
main_analysis_data<-read.csv("./Data/full_data_set_11_29_21_unintentional.csv")
################################## set up data set ################################
#add the intervention dates and time period data
main_analysis_data$Intervention_First_Date<-as.Date(main_analysis_data$Intervention_First_Date)
main_analysis_data$Time_Period_Start<-as.Date(main_analysis_data$Time_Period_Start)
names(main_analysis_data)[which(colnames(main_analysis_data) == "sum_deaths")] <- "imputed_deaths"
################################## set up the Regions ##############################
#set up the regions according to Census: https://www.census.gov/geographies/reference-maps/2010/geo/2010-census-regions-and-divisions-of-the-united-states.html
NE.name <- c("Connecticut","Maine","Massachusetts","New Hampshire",
"Rhode Island","Vermont","New Jersey","New York",
"Pennsylvania")
MW.name <- c("Indiana","Illinois","Michigan","Ohio","Wisconsin",
"Iowa","Kansas","Minnesota","Missouri","Nebraska",
"North Dakota","South Dakota")
S.name <- c("Delaware","District of Columbia","Florida","Georgia",
"Maryland","North Carolina","South Carolina","Virginia",
"West Virginia","Alabama","Kentucky","Mississippi",
"Tennessee","Arkansas","Louisiana","Oklahoma","Texas")
W.name <- c("Arizona","Colorado","Idaho","New Mexico","Montana",
"Utah","Nevada","Wyoming","Alaska","California",
"Hawaii","Oregon","Washington")
region.list <- list(
Northeast=NE.name,
Midwest=MW.name,
South=S.name,
West=W.name)
#initialize vector with "West" and then impute the other regions for the states
main_analysis_data$Region<-rep("West", nrow(main_analysis_data))
for(state in unique(main_analysis_data$State)){
if(state %in% region.list$Northeast){
main_analysis_data$Region[main_analysis_data$State == state]<-"Northeast"
}else if(state %in% region.list$Midwest){
main_analysis_data$Region[main_analysis_data$State == state]<-"Midwest"
}else if(state %in% region.list$South){
main_analysis_data$Region[main_analysis_data$State == state]<-"South"
}
}
#here, we estimate the variance-covariance matrix through the sandwich estimator
#we create a function so that we don't have to keep writing the code:
compute_sd_and_CI_linear_link <- function(Z, log_observed_y, theta, z_value = 1.96, d,
return_full_cov = FALSE){
#this function computes the sandwich estimator for models with linear link function
#Z: the data such that rows are state-time combinations and columns are the different policy measures
#log_observed_y: log of the observed y from the data
#theta: the coefficient values that need to be in order of the columns of cov_data
#z_value the Z-value that corresponds to the CI. We default to 95% CI so we default to 1.96
#d: the number of parameters for a bias correction
#return_full_cov: a TRUE/FALSE boolean that indicates whether to return the full variance-covariance matrix
#compute the middle term of the Sandwich Estimator = sum_{s,t} Z_{s,t}*Z_{s,t}^T * (log(y_{s,t}) - Z_{s,t}^T*theta)^2
middle_term <- matrix(0, nrow = ncol(Z), ncol = ncol(Z))
for(i in 1:nrow(Z)){
#tcrossprod(X) = X*X^T
middle_term <- middle_term + tcrossprod(as.matrix(Z[i,]))*
as.numeric((log_observed_y[i] - t(as.matrix(Z[i,]))%*%theta)^2)
}
#compute C = sum_{s,t} Z_{s,t} * Z_{s,t}^T which is equal to Z^T * Z using matrix multiplication
#we compute C^{-1} by using solve()
C_inverse <- solve(crossprod(Z))
#bias_term = (N)/(N-d) where N = nrow(Z)
bias_term <- nrow(Z)/(nrow(Z) - d)
Sigma_hat <- C_inverse%*%(middle_term)%*%t(C_inverse)*bias_term
#we obtain the standard errors of the coefficients by taking the square root of the diagonal of the variance-covariance matrix.
se_of_coefficients <- sqrt(diag(Sigma_hat))
#find the CI for the coefficients
lb_coef <- theta - z_value*(se_of_coefficients)
ub_coef <- theta + z_value*(se_of_coefficients)
#return_data_set returns the:
#lb_coef: lower bounds of the coefficients
#theta: coefficients
#ub_coef: upper bounds of the coefficients
#sd_coef = standard error of coefficients
return_data_set <- data.frame(lb_coef, theta, ub_coef, se_coef = se_of_coefficients)
#if user wants the full variance-covariance matrix, we return both the return_data_set and Sigma_hat
if(return_full_cov){
return(list(return_data_set = return_data_set, var_cov = Sigma_hat))
}else{
return(return_data_set)
}
}
compute_sd_and_CI_logistic_link <- function(Z, observed_y, theta, z_value = 1.96, d,
return_full_cov = FALSE){
#this function computes the sandwich estimator for models with logistic link function
#Z: the data such that rows are state-time combinations and columns are the different policy measures
#log_observed_y: log of the observed y from the data
#theta: the coefficient values that need to be in order of the columns of cov_data
#z_value the Z-value that corresponds to the CI. We default to 95% CI so we default to 1.96
#d: the number of parameters for a bias correction
#return_full_cov: a TRUE/FALSE boolean that indicates whether to return the full variance-covariance matrix
#initialize the C_term and middle_term
C_term <- middle_term <- matrix(0, nrow = ncol(Z), ncol = ncol(Z))
for(i in 1:nrow(Z)){
#p_{s,t} = expit(z_{s,t}^T Theta)
p_st <- expit(t(as.matrix(Z[i,]))%*%theta)
#middle term = sum_{s,t} Z_{s,t} * Z_{s,t}^T * (Y_{s,t} - p_{s,t})^2
middle_term <- middle_term + tcrossprod(as.matrix(Z[i,]))*
as.numeric((observed_y[i] - p_st)^2)
#C_term = sum_{s,t} Z_{s,t} * Z_{s,t}^T * (p_{s,t}*(1-p_{s,t}))
C_term <- C_term + tcrossprod(as.matrix(Z[i,]))*as.numeric(p_st*(1 - p_st))
}
C_inverse <- solve(C_term)
#bias_term = (N)/(N-d) where N = nrow(Z)
bias_term <- nrow(Z)/(nrow(Z) - d)
Sigma_hat <- C_inverse%*%as.matrix(middle_term)%*%t(C_inverse)*bias_term
#we obtain the standard errors by taking the square root of the diagonal of the variance-covariance matrix.
se_of_coefficients <- sqrt(diag(Sigma_hat))
#find the CI for the coefficients
lb_coef <- theta - z_value*(se_of_coefficients)
ub_coef <- theta + z_value*(se_of_coefficients)
#return_data_set returns the:
#lb_coef: lower bounds of the coefficients
#theta: coefficients
#ub_coef: upper bounds of the coefficients
#exp_lb: exponential of lower bounds of the coefficients
#exp_coef = exponential of theta
#exp_ub: exponential of upper bounds of the coefficients
#sd_coef = standard error of coefficients
return_data_set <- data.frame(lb_coef,
theta,
ub_coef,
exp_lb = exp(lb_coef),
exp_coef = exp(theta),
exp_ub = exp(ub_coef),
se_coef = se_of_coefficients)
#if user wants the full variance-covariance matrix, we return both the return_data_set and Sigma_hat
if(return_full_cov){
return(list(return_data_set = return_data_set, var_cov = Sigma_hat))
}else{
return(return_data_set)
}
}
attr_death_compute <- function(data, coef_data, tx_name = NULL,
var_cov_mat_beta = NULL, prob_of_od_var = "prop_dead"){
attr_table <- data.frame(matrix(NA, nrow = unique(data$Time_Period_ID), ncol = 4))
#filter data so that it's only states where there was treatment
treated_data <- data %>%
filter(Intervention_Redefined > 0)
for(time in unique(treated_data$Time_Period_ID)){
#filter treated_data to time period t
time_data <- treated_data %>%
filter(Time_Period_ID == time)
#obtain the population
pop <- time_data$population
#obtain the estimated probability had intervention not occurred
#here, we compute x^T*beta where x is a vector of the covariates and beta is the corresponding coefficients
est_prob_no_int <- time_data[,prob_of_od_var]*exp(- as.matrix(time_data[,tx_name])%*%as.matrix(coef_data[tx_name, "estimate"]))
#estimated number of OD had intervention not occurred
n_od_no_int <- pop*est_prob_no_int
#obtain LB and UB
if(length(tx_name) == 1){
#if there is only one variable corresponding to treatment, we don't need to account for covariances
est_prob_no_int_lb <- time_data[,prob_of_od_var]*exp(- as.matrix(time_data[,tx_name])%*%as.matrix(coef_data[tx_name, "conf.low"]))
n_attr_od_lb <- sum(pop*est_prob_no_int_lb) - sum(time_data$imputed_deaths)
est_prob_no_int_ub <- time_data[,prob_of_od_var]*exp(-as.matrix(time_data[,tx_name])%*%as.matrix(coef_data[tx_name, "conf.high"]))
n_attr_od_ub <- sum(pop*est_prob_no_int_ub) - sum(time_data$imputed_deaths)
}else{
#otherwise, if there are more than 1 treatment parameters, we need to account for covariances
#here we assume its two treatment parameters
jacobian_first_entry <- sum(time_data$population*time_data[,prob_of_od_var]*
exp(-as.matrix(time_data[,tx_name])%*%as.matrix(coef_data[tx_name,"estimate"]))*
time_data[,tx_name[1]])
jacobian_second_entry <- sum(time_data$population*time_data[,prob_of_od_var]*
exp(-as.matrix(time_data[,tx_name])%*%as.matrix(coef_data[tx_name,"estimate"]))*
time_data[,tx_name[2]])
jacobian_mat <- matrix(c(jacobian_first_entry, jacobian_second_entry), nrow = 1)
computed_sd <- sqrt(tcrossprod(jacobian_mat%*%var_cov_mat_beta, jacobian_mat))
n_attr_od_lb <- sum(n_od_no_int - time_data$imputed_deaths) - 1.96*computed_sd
n_attr_od_ub <- (sum(n_od_no_int - time_data$imputed_deaths)) + 1.96*computed_sd
}
attr_table[time,] <- c(time, sum(n_od_no_int - time_data$imputed_deaths),
(n_attr_od_lb) ,
(n_attr_od_ub) )
}
colnames(attr_table) <- c("Time_Period", "attr_deaths", "attr_deaths_lb", "attr_deaths_ub")
attr_table
}
main_analysis_data$prop_dead <- main_analysis_data$imputed_deaths/main_analysis_data$population
#create the dataset for the event study to check for pre-trend analysis
#First create dataset that indicates the intervention time of each state.
#We then indicate the ID of the six-month time period that the intervention time is in so that we can use this to compare time points
time_data_int <- main_analysis_data %>%
#group by the state
group_by(State) %>%
#find the time interval ID for the intervention time
#we do this by finding the six-month time ID of the intervention date using floor_date(, "6 months")
summarise(intervention_time_id = ifelse(floor_date(Intervention_First_Date, "6 months") == Time_Period_Start, Time_Period_ID, NA)) %>%
#filter out the other six-month time periods that aren't the intervention date
filter(!is.na(intervention_time_id))
#merge the time_data_int with the main dataset by state
merged_main_time_data_int <- merge(main_analysis_data, time_data_int, by = "State", all.x = TRUE)
#create the columns that associate with the periods before the intervention
#we call these negative periods since it is negative relative to the intervention time
#neg_periods_df is a dataset of indicators of x periods before intervention, where rows are equal to number of periods before intervention
#the maximum number of periods pre-intervention is equal to max(intervention_time_id) - 1
neg_periods_df <- sapply(1:(max(merged_main_time_data_int$intervention_time_id, na.rm = TRUE) - 1),
#the indicator for x periods before intervention is equal to 1 if the time ID of intervention minus time ID is equal to x
#otherwise, it is 0
function(x){ifelse(merged_main_time_data_int$intervention_time_id -
merged_main_time_data_int$Time_Period_ID == x, 1, 0)})
#create the column names for neg_periods_df of the form: "neg_x_pd"
colnames(neg_periods_df) <- sapply(1:(max(merged_main_time_data_int$intervention_time_id, na.rm = TRUE) - 1),
function(x){paste("neg_", x, "_pd", sep = "")})
#add in the state and time ID columns to neg_periods_df
neg_periods_df <- cbind(neg_periods_df, "State" = merged_main_time_data_int$State,
"Time_Period_ID" = merged_main_time_data_int$Time_Period_ID)
#for Hawaii, impute a 0 because it is NA right now
neg_periods_df[neg_periods_df[,"State"] == "Hawaii", 1:34] <- 0
#create the columns that associate with the periods after the intervention
#we call these positive periods since it is positive relative to the intervention time
#pos_periods_df is a dataset of indicators of x periods after intervention, where rows are equal to number of periods after intervention
#the max number of periods after the intervention is determined by the maximum Time ID minus the earliest time period of the intervention
#the period 0 is associated with intervention time
pos_periods_df <- sapply(0:(max(merged_main_time_data_int$Time_Period_ID) -
min(merged_main_time_data_int$intervention_time_id, na.rm = TRUE)),
function(x){ifelse(merged_main_time_data_int$Time_Period_ID -
merged_main_time_data_int$intervention_time_id == x, 1, 0)})
#create the column names for pos_periods_df of the form: "pos_x_pd"
colnames(pos_periods_df) <- sapply(0:(max(merged_main_time_data_int$Time_Period_ID) -
min(merged_main_time_data_int$intervention_time_id, na.rm = TRUE)),
function(x){paste("pos_", x, "_pd", sep = "")})
#add in the state and time ID columns
pos_periods_df <- cbind(pos_periods_df, "State" = merged_main_time_data_int$State,
"Time_Period_ID" = merged_main_time_data_int$Time_Period_ID)
#for Hawaii, impute a 0 because it is NA right now
pos_periods_df[pos_periods_df[,"State"] == "Hawaii", 1:40] <- 0
#merge the columns of indicators for before and after the intervention with the main analysis data to create the dataset for event study
sensitivity_anlys_event_study_data <- merge(main_analysis_data,
neg_periods_df, by = c("State", "Time_Period_ID"))
sensitivity_anlys_event_study_data <- merge(sensitivity_anlys_event_study_data,
pos_periods_df, by = c("State", "Time_Period_ID"))
#change the indicator values to numeric type
#we first get the index of the new columns
neg_1_index <- which(colnames(sensitivity_anlys_event_study_data) == "neg_1_pd")
pos_39_index <- which(colnames(sensitivity_anlys_event_study_data) == "pos_39_pd")
#then turn the values to numeric
sensitivity_anlys_event_study_data[, neg_1_index:pos_39_index] <- apply(sensitivity_anlys_event_study_data[, neg_1_index:pos_39_index],
2, as.numeric)
#compute the proportion of people who died from drug overdose
main_analysis_data$prop_dead <- main_analysis_data$imputed_deaths/main_analysis_data$population
#add a variable as a measure of number of states with DIH prosecution
main_analysis_data <- main_analysis_data %>%
group_by(Time_Period_Start) %>%
mutate(num_states_w_intervention = sum(Intervention_Redefined))
#fit an OLS with smoothed time effects
main_analysis_model_log_smoothed_time<-gam(log(prop_dead)~ State +
s(Time_Period_ID, bs = "cr", by = as.factor(Region)) +
Naloxone_Pharmacy_Yes_Redefined +
Naloxone_Pharmacy_No_Redefined +
Medical_Marijuana_Redefined +
Recreational_Marijuana_Redefined +
GSL_Redefined +
PDMP_Redefined +
Medicaid_Expansion_Redefined +
Intervention_Redefined +
num_states_w_intervention,
data = main_analysis_data)
summary(main_analysis_model_log_smoothed_time)
##
## Family: gaussian
## Link function: identity
##
## Formula:
## log(prop_dead) ~ State + s(Time_Period_ID, bs = "cr", by = as.factor(Region)) +
## Naloxone_Pharmacy_Yes_Redefined + Naloxone_Pharmacy_No_Redefined +
## Medical_Marijuana_Redefined + Recreational_Marijuana_Redefined +
## GSL_Redefined + PDMP_Redefined + Medicaid_Expansion_Redefined +
## Intervention_Redefined + num_states_w_intervention
##
## Parametric coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -9.5901401 0.2059517 -46.565 < 2e-16 ***
## StateAlaska 0.1642934 0.0743934 2.208 0.027331 *
## StateArizona 0.2883356 0.0677554 4.256 2.19e-05 ***
## StateArkansas -0.4868587 0.0667523 -7.294 4.40e-13 ***
## StateCalifornia -0.1561194 0.0743859 -2.099 0.035967 *
## StateColorado 0.0361092 0.0741355 0.487 0.626265
## StateConnecticut 0.2301798 0.0709034 3.246 0.001189 **
## StateDelaware 0.2169779 0.0680081 3.190 0.001443 **
## StateFlorida 0.3085987 0.0668104 4.619 4.11e-06 ***
## StateGeorgia 0.0002845 0.0667934 0.004 0.996602
## StateHawaii -0.3783288 0.0729745 -5.184 2.40e-07 ***
## StateIdaho -0.1242470 0.0668007 -1.860 0.063043 .
## StateIllinois 0.1612842 0.0678109 2.378 0.017483 *
## StateIndiana 0.0819325 0.0663673 1.235 0.217156
## StateIowa -0.6769694 0.0665778 -10.168 < 2e-16 ***
## StateKansas -0.2307812 0.0660706 -3.493 0.000489 ***
## StateKentucky 0.7028805 0.0667944 10.523 < 2e-16 ***
## StateLouisiana 0.3407484 0.0659833 5.164 2.67e-07 ***
## StateMaine 0.0880973 0.0738853 1.192 0.233271
## StateMaryland -1.5131006 0.0677976 -22.318 < 2e-16 ***
## StateMassachusetts -0.0833781 0.0673507 -1.238 0.215879
## StateMichigan 0.0099350 0.0683957 0.145 0.884522
## StateMinnesota -0.6430300 0.0699022 -9.199 < 2e-16 ***
## StateMississippi -0.0490213 0.0660266 -0.742 0.457907
## StateMissouri 0.1675455 0.0682763 2.454 0.014219 *
## StateMontana -0.4401972 0.0702845 -6.263 4.65e-10 ***
## StateNebraska -0.8742551 0.0670928 -13.031 < 2e-16 ***
## StateNevada 0.4597148 0.0717591 6.406 1.87e-10 ***
## StateNew Hampshire 0.1476303 0.0671583 2.198 0.028051 *
## StateNew Jersey 0.0807155 0.0679571 1.188 0.235082
## StateNew Mexico 0.6453680 0.0729173 8.851 < 2e-16 ***
## StateNew York -0.1394535 0.0683701 -2.040 0.041518 *
## StateNorth Carolina 0.2339425 0.0658631 3.552 0.000392 ***
## StateNorth Dakota -1.1371560 0.0664632 -17.110 < 2e-16 ***
## StateOhio 0.4525904 0.0669815 6.757 1.86e-11 ***
## StateOklahoma 0.4639698 0.0664674 6.980 4.04e-12 ***
## StateOregon -0.2715116 0.0738537 -3.676 0.000243 ***
## StatePennsylvania 0.5608113 0.0668628 8.387 < 2e-16 ***
## StateRhode Island -0.2786824 0.0693850 -4.016 6.14e-05 ***
## StateSouth Carolina 0.2077684 0.0664045 3.129 0.001781 **
## StateSouth Dakota -1.0271060 0.0667936 -15.377 < 2e-16 ***
## StateTennessee 0.4686941 0.0656809 7.136 1.36e-12 ***
## StateTexas -0.0198795 0.0667433 -0.298 0.765850
## StateUtah -0.0957968 0.0660724 -1.450 0.147257
## StateVermont -0.1694961 0.0696873 -2.432 0.015097 *
## StateVirginia -0.0323098 0.0659867 -0.490 0.624444
## StateWashington 0.0451432 0.0752379 0.600 0.548573
## StateWest Virginia 0.7907588 0.0668153 11.835 < 2e-16 ***
## StateWisconsin 0.0062720 0.0661277 0.095 0.924446
## StateWyoming -0.0218261 0.0660516 -0.330 0.741103
## Naloxone_Pharmacy_Yes_Redefined -0.0802046 0.0424748 -1.888 0.059138 .
## Naloxone_Pharmacy_No_Redefined -0.0014359 0.0385076 -0.037 0.970259
## Medical_Marijuana_Redefined 0.1928104 0.0306582 6.289 3.95e-10 ***
## Recreational_Marijuana_Redefined -0.1113138 0.0487806 -2.282 0.022603 *
## GSL_Redefined 0.0536153 0.0315983 1.697 0.089900 .
## PDMP_Redefined -0.1526621 0.0246836 -6.185 7.58e-10 ***
## Medicaid_Expansion_Redefined 0.0909280 0.0301320 3.018 0.002581 **
## Intervention_Redefined -0.0253718 0.0244012 -1.040 0.298574
## num_states_w_intervention -0.0042330 0.0069748 -0.607 0.543993
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Approximate significance of smooth terms:
## edf Ref.df F p-value
## s(Time_Period_ID):as.factor(Region)Midwest 4.678 5.750 12.35 <2e-16 ***
## s(Time_Period_ID):as.factor(Region)Northeast 8.454 8.914 17.38 <2e-16 ***
## s(Time_Period_ID):as.factor(Region)South 6.360 7.506 12.53 <2e-16 ***
## s(Time_Period_ID):as.factor(Region)West 3.027 3.846 10.72 <2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## R-sq.(adj) = 0.84 Deviance explained = 84.6%
## GCV = 0.089657 Scale est. = 0.086002 n = 2000
gam.check(main_analysis_model_log_smoothed_time, page = 1)
##
## Method: GCV Optimizer: magic
## Smoothing parameter selection converged after 5 iterations.
## The RMS GCV score gradient at convergence was 1.113005e-06 .
## The Hessian was positive definite.
## Model rank = 95 / 95
##
## Basis dimension (k) checking results. Low p-value (k-index<1) may
## indicate that k is too low, especially if edf is close to k'.
##
## k' edf k-index p-value
## s(Time_Period_ID):as.factor(Region)Midwest 9.00 4.68 1.05 0.99
## s(Time_Period_ID):as.factor(Region)Northeast 9.00 8.45 1.05 0.99
## s(Time_Period_ID):as.factor(Region)South 9.00 6.36 1.05 0.99
## s(Time_Period_ID):as.factor(Region)West 9.00 3.03 1.05 0.99
#examine fitted values
summary(fitted(main_analysis_model_log_smoothed_time))
## Min. 1st Qu. Median Mean 3rd Qu. Max.
## -12.316 -10.207 -9.739 -9.796 -9.343 -8.122
hist(fitted(main_analysis_model_log_smoothed_time))
#smoothed effects
plot(main_analysis_model_log_smoothed_time, pages = 1)
#compute the full dataset including basis functions (these are the functions used to compute smoothed time effects)
full_df_w_basis_functions_log_smoothed_time <- as.matrix(data.frame(predict(main_analysis_model_log_smoothed_time, type = "lpmatrix")))
#estimate the 95% CI and SE
coefficient_values_log_smoothed_time <- coef(main_analysis_model_log_smoothed_time)
#obtain the table of coefficients and standard errors
main_analysis_sd_and_ci_log_smoothed_time <- compute_sd_and_CI_linear_link(full_df_w_basis_functions_log_smoothed_time,
log(main_analysis_data$prop_dead),
coefficient_values_log_smoothed_time,
d = ncol(full_df_w_basis_functions_log_smoothed_time))
colnames(main_analysis_sd_and_ci_log_smoothed_time) <- c("conf.low", "estimate", "conf.high", "sd")
#create a table for the coefficients
main_analysis_sd_and_ci_log_smoothed_time$term <- rownames(main_analysis_sd_and_ci_log_smoothed_time)
#create a plot for the coefficients
main_analysis_sd_and_ci_log_smoothed_time$ci_95 <-
paste("95% CI = (", format(round(main_analysis_sd_and_ci_log_smoothed_time$conf.low, 3), nsmall = 3), ", ",
format(round(main_analysis_sd_and_ci_log_smoothed_time$conf.high, 3), nsmall = 3), ")", sep = "")
dwplot(main_analysis_sd_and_ci_log_smoothed_time[51:59,], colour = "black") +
theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank(),
panel.background = element_blank(), axis.line = element_line(colour = "black")) +
geom_vline(aes(xintercept = 0), linetype = "dashed") +
labs(y = "Term", x = "Coefficients and 95% Confidence Intervals",
title = "Coefficient of Analysis With Smoothed Time Effects,
Exposure Intervention") +
scale_color_grey() +
geom_text(main_analysis_sd_and_ci_log_smoothed_time[51:59,],
mapping = aes(label = format(round(estimate, 3), nsmall = 3), x = 0.55, y = 9:1), size = 3) +
geom_text(main_analysis_sd_and_ci_log_smoothed_time[51:59,],
mapping = aes(label = ci_95, x = 0.9, y = 9:1), size = 3) +
xlim(-.3, 1.1)
#create the table of risk ratios of the coefficients
main_analysis_coef_of_interest_RR <- main_analysis_sd_and_ci_log_smoothed_time[51:59,]
main_analysis_coef_of_interest_RR$est_RR <- round(exp(main_analysis_coef_of_interest_RR$estimate), 3)
main_analysis_coef_of_interest_RR$est_RR_lb <- round(exp(main_analysis_coef_of_interest_RR$conf.low), 3)
main_analysis_coef_of_interest_RR$est_RR_ub <- round(exp(main_analysis_coef_of_interest_RR$conf.high), 3)
main_analysis_coef_of_interest_RR[,c("est_RR", "est_RR_lb", "est_RR_ub")]
## est_RR est_RR_lb est_RR_ub
## Naloxone_Pharmacy_Yes_Redefined 0.923 0.863 0.987
## Naloxone_Pharmacy_No_Redefined 0.999 0.936 1.065
## Medical_Marijuana_Redefined 1.213 1.133 1.298
## Recreational_Marijuana_Redefined 0.895 0.832 0.962
## GSL_Redefined 1.055 0.999 1.115
## PDMP_Redefined 0.858 0.813 0.906
## Medicaid_Expansion_Redefined 1.095 1.040 1.153
## Intervention_Redefined 0.975 0.931 1.021
## num_states_w_intervention 0.996 0.979 1.013
date_data <- main_analysis_data[, c("Time_Period_ID", "Time_Period_Start")]
date_data <- date_data[!duplicated(date_data),]
attr_deaths_est_log_smoothed_time <- attr_death_compute(main_analysis_data, main_analysis_sd_and_ci_log_smoothed_time,
tx_name = "Intervention_Redefined")
attr_deaths_est_log_smoothed_time <- merge(attr_deaths_est_log_smoothed_time, date_data, by.x = "Time_Period", by.y = "Time_Period_ID")
ggplot(attr_deaths_est_log_smoothed_time, aes(x = Time_Period_Start)) +
# geom_point(aes(y = attr_deaths)) +
geom_line(aes(y = attr_deaths, linetype = "Estimate")) +
# geom_point(aes(y = attr_deaths_lb)) +
geom_line(aes(y = attr_deaths_lb, linetype = "95% CI")) +
# geom_point(aes(y = attr_deaths_ub)) +
geom_line(aes(y = attr_deaths_ub, linetype = "95% CI")) +
labs(x = "Date", y = "Lives Saved",
title = "Estimated Number of Lives Saved Using Smoothed Time Effects,
Log Probability of Drug Overdose Death",
linetype = "") +
theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank(),
panel.background = element_blank(), axis.line = element_line(colour = "black")) +
scale_linetype_manual(values = c("dashed", "solid"))
#create a formula for the gam model which includes the state effects, smoothed time effects, policy measures,
#the periods before the intervention (excluding 1 period and 34 periods before intervention)
#the intervention period, and the periods after the intervention
formula_event_study_log_smoothed_time <- formula(paste("log(prop_dead) ~ State +
s(Time_Period_ID, bs = 'cr', by = as.factor(Region)) +
Naloxone_Pharmacy_Yes_Redefined +
Naloxone_Pharmacy_No_Redefined +
Medical_Marijuana_Redefined +
Recreational_Marijuana_Redefined +
GSL_Redefined +
PDMP_Redefined +
Medicaid_Expansion_Redefined +",
paste(sapply(2:(max(merged_main_time_data_int$intervention_time_id, na.rm = TRUE)-2),
function(x)paste("neg_", x, "_pd", sep = "")), collapse = "+"),
"+",
paste(sapply(0:(max(merged_main_time_data_int$Time_Period_ID) -
min(merged_main_time_data_int$intervention_time_id, na.rm = TRUE)),
function(x)paste("pos_", x, "_pd", sep = "")), collapse = "+")))
#run the gam model
sensitivity_anlys_event_study_model_log_smoothed_time<-gam(formula_event_study_log_smoothed_time,
data = sensitivity_anlys_event_study_data)
# summary(sensitivity_anlys_event_study_model_log_smoothed_time)
#compute the full dataset including basis functions
full_df_w_basis_functions_sensitivity_anlys_event_study_log_smoothed_time <-
data.frame(predict(sensitivity_anlys_event_study_model_log_smoothed_time, type = "lpmatrix"))
#estimate the 95% CI and SE
coefficient_values_sensitivity_anlys_event_study_log_smoothed_time <- coef(sensitivity_anlys_event_study_model_log_smoothed_time)
sensitivity_anlys_event_study_sd_and_ci_log_smoothed_time <-
compute_sd_and_CI_linear_link(as.matrix(full_df_w_basis_functions_sensitivity_anlys_event_study_log_smoothed_time),
log(sensitivity_anlys_event_study_data$prop_dead),
coefficient_values_sensitivity_anlys_event_study_log_smoothed_time,
d = ncol(full_df_w_basis_functions_sensitivity_anlys_event_study_log_smoothed_time))
(sensitivity_anlys_event_study_sd_and_ci_log_smoothed_time)
## lb_coef theta
## (Intercept) -9.730462237 -9.607755377
## StateAlaska -0.049423874 0.101715259
## StateArizona 0.144496017 0.257887425
## StateArkansas -0.626540747 -0.516038962
## StateCalifornia -0.197605999 -0.062450464
## StateColorado -0.121794519 0.002447472
## StateConnecticut 0.087500824 0.222377470
## StateDelaware 0.018646885 0.155066167
## StateFlorida 0.292773480 0.412924559
## StateGeorgia 0.012368305 0.127609736
## StateHawaii -0.610276068 -0.467150486
## StateIdaho -0.276925613 -0.159066931
## StateIllinois 0.103069959 0.233581709
## StateIndiana -0.024961041 0.076071921
## StateIowa -0.734431274 -0.624312899
## StateKansas -0.322260429 -0.224552743
## StateKentucky 0.591020827 0.686996193
## StateLouisiana 0.292097065 0.394888461
## StateMaine -0.065746614 0.074857829
## StateMaryland -1.675419104 -1.451851831
## StateMassachusetts -0.302002541 -0.083062578
## StateMichigan -0.046464767 0.050576690
## StateMinnesota -0.754905369 -0.645251331
## StateMississippi -0.210742379 -0.087268616
## StateMissouri 0.085057722 0.185702542
## StateMontana -0.510770699 -0.394331541
## StateNebraska -1.018596137 -0.912234589
## StateNevada 0.384505286 0.495672978
## StateNew Hampshire 0.010674384 0.124591896
## StateNew Jersey -0.010708998 0.145727757
## StateNew Mexico 0.486968538 0.622210509
## StateNew York -0.252658248 -0.132279799
## StateNorth Carolina 0.186121130 0.274412387
## StateNorth Dakota -1.343680088 -1.172707883
## StateOhio 0.453871910 0.589151698
## StateOklahoma 0.336286341 0.453117480
## StateOregon -0.416865124 -0.299059866
## StatePennsylvania 0.559441441 0.677459526
## StateRhode Island -0.589873781 -0.312527565
## StateSouth Carolina 0.045203737 0.161134027
## StateSouth Dakota -1.218317372 -1.084870928
## StateTennessee 0.382305713 0.470729691
## StateTexas -0.033873294 0.087487525
## StateUtah -0.241196453 -0.054041819
## StateVermont -0.335254127 -0.191087648
## StateVirginia -0.069301461 0.031405848
## StateWashington -0.096970886 0.024086450
## StateWest Virginia 0.616272337 0.763862290
## StateWisconsin -0.041315108 0.055350995
## StateWyoming -0.158016158 -0.026862195
## Naloxone_Pharmacy_Yes_Redefined -0.126819690 -0.061319716
## Naloxone_Pharmacy_No_Redefined -0.061827052 0.003373398
## Medical_Marijuana_Redefined 0.129050881 0.199070518
## Recreational_Marijuana_Redefined -0.175467189 -0.098626529
## GSL_Redefined 0.007698979 0.063283495
## PDMP_Redefined -0.234541848 -0.178713621
## Medicaid_Expansion_Redefined 0.037370623 0.089195898
## neg_2_pd -0.092349592 0.020524738
## neg_3_pd -0.094755826 0.016757110
## neg_4_pd -0.131739932 -0.010783798
## neg_5_pd -0.130295904 -0.013004948
## neg_6_pd -0.117049085 -0.005778559
## neg_7_pd -0.160704646 -0.044021509
## neg_8_pd -0.245224820 -0.100544147
## neg_9_pd -0.182129446 -0.025158944
## neg_10_pd -0.118922597 0.024031356
## neg_11_pd -0.125372091 0.025412131
## neg_12_pd -0.050114192 0.115435529
## neg_13_pd -0.152868158 0.015177292
## neg_14_pd -0.174761759 -0.015729613
## neg_15_pd -0.248669745 -0.059495317
## neg_16_pd -0.217482729 -0.043767301
## neg_17_pd -0.213246323 -0.060151411
## neg_18_pd -0.196116002 -0.051908104
## neg_19_pd -0.386274566 -0.154364075
## neg_20_pd -0.449689904 -0.193711115
## neg_21_pd -0.299795811 -0.117692721
## neg_22_pd -0.287517173 -0.096423617
## neg_23_pd -0.402605591 -0.116590871
## neg_24_pd -0.522485297 -0.180864324
## neg_25_pd -0.325051268 -0.038788957
## neg_26_pd -0.306147219 0.044238456
## neg_27_pd -0.714236416 -0.213628304
## neg_28_pd -0.669383634 -0.201311333
## neg_29_pd -0.245517266 0.024088229
## neg_30_pd -0.307566464 -0.052096370
## neg_31_pd -0.277751718 -0.054179201
## neg_32_pd -0.179849349 0.032863102
## neg_33_pd -0.216537378 0.109007659
## pos_0_pd -0.119247683 -0.011149669
## pos_1_pd -0.155493502 -0.040113963
## pos_2_pd -0.110679260 -0.005207113
## pos_3_pd -0.149392401 -0.041493922
## pos_4_pd -0.155193824 -0.047126354
## pos_5_pd -0.201852471 -0.092045551
## pos_6_pd -0.211676132 -0.097121035
## pos_7_pd -0.214556474 -0.094523821
## pos_8_pd -0.254768298 -0.140578750
## pos_9_pd -0.281759585 -0.162175014
## pos_10_pd -0.306994676 -0.180016549
## pos_11_pd -0.302439201 -0.183783467
## pos_12_pd -0.317775638 -0.193795841
## pos_13_pd -0.402797602 -0.262320233
## pos_14_pd -0.387924786 -0.247335909
## pos_15_pd -0.387322987 -0.250438341
## pos_16_pd -0.409997721 -0.268524051
## pos_17_pd -0.406130067 -0.263061724
## pos_18_pd -0.400316748 -0.253663735
## pos_19_pd -0.385391286 -0.243608362
## pos_20_pd -0.398712925 -0.247582690
## pos_21_pd -0.434854946 -0.274332431
## pos_22_pd -0.424666741 -0.266749984
## pos_23_pd -0.439924967 -0.273338184
## pos_24_pd -0.487458699 -0.310854938
## pos_25_pd -0.475905008 -0.279921074
## pos_26_pd -0.481207255 -0.289983620
## pos_27_pd -0.541287653 -0.338308822
## pos_28_pd -0.530318463 -0.321207072
## pos_29_pd -0.554552376 -0.334408979
## pos_30_pd -0.513141915 -0.285008197
## pos_31_pd -0.496762582 -0.241545861
## pos_32_pd -0.548421385 -0.298666643
## pos_33_pd -0.636276686 -0.326206090
## pos_34_pd -0.633335400 -0.314276284
## pos_35_pd -0.790567646 -0.514644392
## pos_36_pd -0.805117655 -0.541250661
## pos_37_pd -0.781870141 -0.534987709
## pos_38_pd -0.755748754 -0.513830174
## pos_39_pd -0.844780244 -0.534088308
## s(Time_Period_ID):as.factor(Region)Midwest.1 -0.606300486 -0.494737079
## s(Time_Period_ID):as.factor(Region)Midwest.2 -0.313196740 -0.224608521
## s(Time_Period_ID):as.factor(Region)Midwest.3 0.005273923 0.093175900
## s(Time_Period_ID):as.factor(Region)Midwest.4 0.247618782 0.333711394
## s(Time_Period_ID):as.factor(Region)Midwest.5 0.429783397 0.513087394
## s(Time_Period_ID):as.factor(Region)Midwest.6 0.613350902 0.703168971
## s(Time_Period_ID):as.factor(Region)Midwest.7 0.809910433 0.903771029
## s(Time_Period_ID):as.factor(Region)Midwest.8 1.040189148 1.159056321
## s(Time_Period_ID):as.factor(Region)Midwest.9 0.970869551 1.099762557
## s(Time_Period_ID):as.factor(Region)Northeast.1 -0.858979108 -0.570605352
## s(Time_Period_ID):as.factor(Region)Northeast.2 -0.578240648 -0.351951731
## s(Time_Period_ID):as.factor(Region)Northeast.3 0.135544664 0.305946250
## s(Time_Period_ID):as.factor(Region)Northeast.4 0.069629519 0.234789810
## s(Time_Period_ID):as.factor(Region)Northeast.5 0.184484866 0.353117285
## s(Time_Period_ID):as.factor(Region)Northeast.6 0.438958937 0.593323854
## s(Time_Period_ID):as.factor(Region)Northeast.7 0.822062192 0.987559759
## s(Time_Period_ID):as.factor(Region)Northeast.8 1.245895974 1.424383187
## s(Time_Period_ID):as.factor(Region)Northeast.9 1.039238207 1.204455894
## s(Time_Period_ID):as.factor(Region)South.1 -0.486926536 -0.389922135
## s(Time_Period_ID):as.factor(Region)South.2 -0.208412768 -0.131694247
## s(Time_Period_ID):as.factor(Region)South.3 0.057708791 0.143187991
## s(Time_Period_ID):as.factor(Region)South.4 0.244212376 0.320084111
## s(Time_Period_ID):as.factor(Region)South.5 0.395655933 0.464161403
## s(Time_Period_ID):as.factor(Region)South.6 0.523618590 0.604633801
## s(Time_Period_ID):as.factor(Region)South.7 0.707913502 0.807908570
## s(Time_Period_ID):as.factor(Region)South.8 0.957313615 1.095942563
## s(Time_Period_ID):as.factor(Region)South.9 0.824688260 0.991245595
## s(Time_Period_ID):as.factor(Region)West.1 -0.428787038 -0.308386528
## s(Time_Period_ID):as.factor(Region)West.2 -0.225040794 -0.131400014
## s(Time_Period_ID):as.factor(Region)West.3 -0.021682692 0.069973226
## s(Time_Period_ID):as.factor(Region)West.4 0.154473915 0.230646763
## s(Time_Period_ID):as.factor(Region)West.5 0.275587256 0.357605994
## s(Time_Period_ID):as.factor(Region)West.6 0.374227506 0.468020533
## s(Time_Period_ID):as.factor(Region)West.7 0.467281938 0.566258884
## s(Time_Period_ID):as.factor(Region)West.8 0.602067415 0.719563215
## s(Time_Period_ID):as.factor(Region)West.9 0.585285435 0.704585511
## ub_coef se_coef
## (Intercept) -9.485048518 0.06260554
## StateAlaska 0.252854392 0.07711180
## StateArizona 0.371278832 0.05785276
## StateArkansas -0.405537177 0.05637846
## StateCalifornia 0.072705071 0.06895691
## StateColorado 0.126689462 0.06338877
## StateConnecticut 0.357254117 0.06881462
## StateDelaware 0.291485448 0.06960167
## StateFlorida 0.533075637 0.06130157
## StateGeorgia 0.242851168 0.05879665
## StateHawaii -0.324024904 0.07302326
## StateIdaho -0.041208248 0.06013198
## StateIllinois 0.364093459 0.06658763
## StateIndiana 0.177104884 0.05154743
## StateIowa -0.514194525 0.05618284
## StateKansas -0.126845056 0.04985086
## StateKentucky 0.782971560 0.04896702
## StateLouisiana 0.497679858 0.05244459
## StateMaine 0.215462271 0.07173696
## StateMaryland -1.228284558 0.11406494
## StateMassachusetts 0.135877385 0.11170406
## StateMichigan 0.147618148 0.04951095
## StateMinnesota -0.535597293 0.05594594
## StateMississippi 0.036205147 0.06299682
## StateMissouri 0.286347363 0.05134940
## StateMontana -0.277892384 0.05940773
## StateNebraska -0.805873041 0.05426610
## StateNevada 0.606840670 0.05671821
## StateNew Hampshire 0.238509408 0.05812118
## StateNew Jersey 0.302164511 0.07981467
## StateNew Mexico 0.757452480 0.06900101
## StateNew York -0.011901350 0.06141758
## StateNorth Carolina 0.362703644 0.04504656
## StateNorth Dakota -1.001735677 0.08723072
## StateOhio 0.724431486 0.06902030
## StateOklahoma 0.569948619 0.05960772
## StateOregon -0.181254607 0.06010472
## StatePennsylvania 0.795477611 0.06021331
## StateRhode Island -0.035181348 0.14150317
## StateSouth Carolina 0.277064316 0.05914811
## StateSouth Dakota -0.951424485 0.06808492
## StateTennessee 0.559153670 0.04511427
## StateTexas 0.208848344 0.06191879
## StateUtah 0.133112815 0.09548706
## StateVermont -0.046921170 0.07355433
## StateVirginia 0.132113157 0.05138128
## StateWashington 0.145143787 0.06176395
## StateWest Virginia 0.911452242 0.07530100
## StateWisconsin 0.152017097 0.04931944
## StateWyoming 0.104291768 0.06691529
## Naloxone_Pharmacy_Yes_Redefined 0.004180259 0.03341835
## Naloxone_Pharmacy_No_Redefined 0.068573849 0.03326554
## Medical_Marijuana_Redefined 0.269090156 0.03572430
## Recreational_Marijuana_Redefined -0.021785870 0.03920442
## GSL_Redefined 0.118868012 0.02835945
## PDMP_Redefined -0.122885394 0.02848379
## Medicaid_Expansion_Redefined 0.141021172 0.02644147
## neg_2_pd 0.133399068 0.05758894
## neg_3_pd 0.128270046 0.05689435
## neg_4_pd 0.110172337 0.06171231
## neg_5_pd 0.104286007 0.05984232
## neg_6_pd 0.105491968 0.05677068
## neg_7_pd 0.072661628 0.05953221
## neg_8_pd 0.044136525 0.07381667
## neg_9_pd 0.131811558 0.08008699
## neg_10_pd 0.166985310 0.07293569
## neg_11_pd 0.176196353 0.07693073
## neg_12_pd 0.280985250 0.08446414
## neg_13_pd 0.183222741 0.08573747
## neg_14_pd 0.143302533 0.08113885
## neg_15_pd 0.129679111 0.09651757
## neg_16_pd 0.129948127 0.08863032
## neg_17_pd 0.092943501 0.07810965
## neg_18_pd 0.092299794 0.07357546
## neg_19_pd 0.077546416 0.11832168
## neg_20_pd 0.062267675 0.13060142
## neg_21_pd 0.064410369 0.09290974
## neg_22_pd 0.094669940 0.09749671
## neg_23_pd 0.169423848 0.14592588
## neg_24_pd 0.160756649 0.17429641
## neg_25_pd 0.247473355 0.14605220
## neg_26_pd 0.394624131 0.17876820
## neg_27_pd 0.286979809 0.25541230
## neg_28_pd 0.266760967 0.23881240
## neg_29_pd 0.293693724 0.13755382
## neg_30_pd 0.203373723 0.13034188
## neg_31_pd 0.169393317 0.11406761
## neg_32_pd 0.245575554 0.10852676
## neg_33_pd 0.434552697 0.16609441
## pos_0_pd 0.096948344 0.05515205
## pos_1_pd 0.075265577 0.05886711
## pos_2_pd 0.100265034 0.05381232
## pos_3_pd 0.066404557 0.05505024
## pos_4_pd 0.060941117 0.05513646
## pos_5_pd 0.017761369 0.05602394
## pos_6_pd 0.017434062 0.05844648
## pos_7_pd 0.025508832 0.06124115
## pos_8_pd -0.026389202 0.05825997
## pos_9_pd -0.042590442 0.06101254
## pos_10_pd -0.053038422 0.06478476
## pos_11_pd -0.065127733 0.06053864
## pos_12_pd -0.069816043 0.06325500
## pos_13_pd -0.121842865 0.07167213
## pos_14_pd -0.106747033 0.07172902
## pos_15_pd -0.113553695 0.06983911
## pos_16_pd -0.127050380 0.07218044
## pos_17_pd -0.119993382 0.07299405
## pos_18_pd -0.107010723 0.07482297
## pos_19_pd -0.101825438 0.07233823
## pos_20_pd -0.096452455 0.07710726
## pos_21_pd -0.113809917 0.08189924
## pos_22_pd -0.108833226 0.08056977
## pos_23_pd -0.106751400 0.08499326
## pos_24_pd -0.134251176 0.09010396
## pos_25_pd -0.083937141 0.09999180
## pos_26_pd -0.098759986 0.09756308
## pos_27_pd -0.135329992 0.10356063
## pos_28_pd -0.112095682 0.10668949
## pos_29_pd -0.114265581 0.11231806
## pos_30_pd -0.056874479 0.11639475
## pos_31_pd 0.013670859 0.13021261
## pos_32_pd -0.048911901 0.12742589
## pos_33_pd -0.016135493 0.15819928
## pos_34_pd 0.004782832 0.16278526
## pos_35_pd -0.238721139 0.14077717
## pos_36_pd -0.277383666 0.13462602
## pos_37_pd -0.288105278 0.12596042
## pos_38_pd -0.271911593 0.12342785
## pos_39_pd -0.223396372 0.15851629
## s(Time_Period_ID):as.factor(Region)Midwest.1 -0.383173672 0.05692011
## s(Time_Period_ID):as.factor(Region)Midwest.2 -0.136020302 0.04519807
## s(Time_Period_ID):as.factor(Region)Midwest.3 0.181077878 0.04484795
## s(Time_Period_ID):as.factor(Region)Midwest.4 0.419804007 0.04392480
## s(Time_Period_ID):as.factor(Region)Midwest.5 0.596391391 0.04250204
## s(Time_Period_ID):as.factor(Region)Midwest.6 0.792987040 0.04582555
## s(Time_Period_ID):as.factor(Region)Midwest.7 0.997631624 0.04788806
## s(Time_Period_ID):as.factor(Region)Midwest.8 1.277923493 0.06064652
## s(Time_Period_ID):as.factor(Region)Midwest.9 1.228655562 0.06576174
## s(Time_Period_ID):as.factor(Region)Northeast.1 -0.282231597 0.14712947
## s(Time_Period_ID):as.factor(Region)Northeast.2 -0.125662815 0.11545353
## s(Time_Period_ID):as.factor(Region)Northeast.3 0.476347836 0.08693958
## s(Time_Period_ID):as.factor(Region)Northeast.4 0.399950100 0.08426545
## s(Time_Period_ID):as.factor(Region)Northeast.5 0.521749703 0.08603695
## s(Time_Period_ID):as.factor(Region)Northeast.6 0.747688772 0.07875761
## s(Time_Period_ID):as.factor(Region)Northeast.7 1.153057326 0.08443753
## s(Time_Period_ID):as.factor(Region)Northeast.8 1.602870399 0.09106490
## s(Time_Period_ID):as.factor(Region)Northeast.9 1.369673580 0.08429474
## s(Time_Period_ID):as.factor(Region)South.1 -0.292917733 0.04949204
## s(Time_Period_ID):as.factor(Region)South.2 -0.054975725 0.03914210
## s(Time_Period_ID):as.factor(Region)South.3 0.228667190 0.04361184
## s(Time_Period_ID):as.factor(Region)South.4 0.395955845 0.03871007
## s(Time_Period_ID):as.factor(Region)South.5 0.532666874 0.03495177
## s(Time_Period_ID):as.factor(Region)South.6 0.685649013 0.04133429
## s(Time_Period_ID):as.factor(Region)South.7 0.907903637 0.05101789
## s(Time_Period_ID):as.factor(Region)South.8 1.234571510 0.07072905
## s(Time_Period_ID):as.factor(Region)South.9 1.157802930 0.08497823
## s(Time_Period_ID):as.factor(Region)West.1 -0.187986018 0.06142883
## s(Time_Period_ID):as.factor(Region)West.2 -0.037759235 0.04777591
## s(Time_Period_ID):as.factor(Region)West.3 0.161629143 0.04676322
## s(Time_Period_ID):as.factor(Region)West.4 0.306819612 0.03886370
## s(Time_Period_ID):as.factor(Region)West.5 0.439624731 0.04184629
## s(Time_Period_ID):as.factor(Region)West.6 0.561813560 0.04785359
## s(Time_Period_ID):as.factor(Region)West.7 0.665235830 0.05049844
## s(Time_Period_ID):as.factor(Region)West.8 0.837059016 0.05994684
## s(Time_Period_ID):as.factor(Region)West.9 0.823885587 0.06086739
# write.csv(format(round(sensitivity_anlys_event_study_sd_and_ci, 3), nsmall = 3), "./Data/event_study_coef_and_ci.csv")
#plot the coefficients for the periods before and after the intervention with 95% CI
#we first have to make a dataset for the plot by compiling the parameter name, coefficient, lower bound, and upper bound
#we also just make the plot with the coefficients corresponding to the periods before and after the intervention
plot_event_study_log_smoothed_time <- sensitivity_anlys_event_study_sd_and_ci_log_smoothed_time %>%
mutate(term = rownames(sensitivity_anlys_event_study_sd_and_ci_log_smoothed_time)) %>%
dplyr::select(term, theta, lb_coef, ub_coef) %>%
filter(term %in% c(sapply(2:(max(merged_main_time_data_int$intervention_time_id, na.rm = TRUE) - 2),
function(x){paste("neg_", x, "_pd", sep = "")}),
sapply(0:(max(merged_main_time_data_int$Time_Period_ID) -
min(merged_main_time_data_int$intervention_time_id, na.rm = TRUE)),
function(x){paste("pos_", x, "_pd", sep = "")})))
#rename the column names to match those needed for the dot and whisker plot
colnames(plot_event_study_log_smoothed_time) <- c("term", "estimate", "conf.low", "conf.high")
#clean up the names of the periods by replacing "pos_" with "+" and "neg_" with "-"
plot_event_study_log_smoothed_time_paper <- plot_event_study_log_smoothed_time
plot_event_study_log_smoothed_time_paper$term <- ifelse(grepl("pos_", plot_event_study_log_smoothed_time_paper$term),
paste("+", extract_numeric(plot_event_study_log_smoothed_time_paper$term),
sep = ""),
paste("-", extract_numeric(plot_event_study_log_smoothed_time_paper$term),
sep = ""))
# pdf("./Figures/pre_tx_trend_ols_model_5_13_22.pdf", width = 5, height = 5)
dwplot(plot_event_study_log_smoothed_time_paper, colour = "black",
vars_order = c(sapply((max(merged_main_time_data_int$Time_Period_ID) -
min(merged_main_time_data_int$intervention_time_id, na.rm = TRUE)):0,
function(x){paste("+", x, sep = "")}),
sapply(2:(max(merged_main_time_data_int$intervention_time_id, na.rm = TRUE) - 2),
function(x){paste("-", x, sep = "")}))) +
theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank(),
panel.background = element_blank(), axis.line = element_line(colour = "black"),
axis.text.x = element_text(size = 5, angle = 45)) +
geom_vline(aes(xintercept = 0), linetype = "dashed") +
labs(y = "Time Relative to Time of Treatment",
x = "Coefficients and 95% Confidence Intervals") +
scale_color_grey() +
coord_flip() +
geom_hline(yintercept = 33, col = "red", linetype = "dashed")
# dev.off()
Now we consider the same analysis but without the pre-treatment periods.
formula_post_tx_log_smoothed_time <- formula(paste("log(prop_dead)~ State +
s(Time_Period_ID, bs = 'cr', by = as.factor(Region)) +
Naloxone_Pharmacy_Yes_Redefined +
Naloxone_Pharmacy_No_Redefined +
Medical_Marijuana_Redefined +
Recreational_Marijuana_Redefined +
GSL_Redefined +
PDMP_Redefined +
Medicaid_Expansion_Redefined +",
paste(sapply(0:(max(merged_main_time_data_int$Time_Period_ID) -
min(merged_main_time_data_int$intervention_time_id, na.rm = TRUE)),
function(x)paste("pos_", x, "_pd", sep = "")), collapse = "+")))
#run the gam model
sensitivity_anlys_post_tx_model_log_smoothed_time<-gam(formula_post_tx_log_smoothed_time,
data = sensitivity_anlys_event_study_data)
# summary(sensitivity_anlys_post_tx_model_log_smoothed_time)
#compute the full dataset including basis functions
full_df_w_basis_functions_sensitivity_anlys_post_tx_log_smoothed_time <-
data.frame(predict(sensitivity_anlys_post_tx_model_log_smoothed_time, type = "lpmatrix"))
#estimate the 95% CI and SD
coefficient_values_sensitivity_anlys_post_tx_log_smoothed_time <- coef(sensitivity_anlys_post_tx_model_log_smoothed_time)
sensitivity_anlys_post_tx_sd_and_ci_log_smoothed_time <-
compute_sd_and_CI_linear_link(as.matrix(full_df_w_basis_functions_sensitivity_anlys_post_tx_log_smoothed_time),
log(sensitivity_anlys_event_study_data$prop_dead),
coefficient_values_sensitivity_anlys_post_tx_log_smoothed_time,
d = ncol(full_df_w_basis_functions_sensitivity_anlys_post_tx_log_smoothed_time))
# sensitivity_anlys_post_tx_sd_and_ci_log_smoothed_time
#plot the coefficients for the periods before and after the intervention with 95% CI
plot_post_tx_log_smoothed_time <- sensitivity_anlys_post_tx_sd_and_ci_log_smoothed_time %>%
mutate(term = rownames(sensitivity_anlys_post_tx_sd_and_ci_log_smoothed_time)) %>%
dplyr::select(term, theta, lb_coef, ub_coef) %>%
filter(term %in% c(sapply(0:(max(merged_main_time_data_int$Time_Period_ID) -
min(merged_main_time_data_int$intervention_time_id, na.rm = TRUE)),
function(x){paste("pos_", x, "_pd", sep = "")})))
#rename the column names to match those needed for the dot and whisker plot
colnames(plot_post_tx_log_smoothed_time) <- c("term", "estimate", "conf.low", "conf.high")
#compute the number of states with x periods post treatment to have an idea of the amount of data
plot_post_tx_log_smoothed_time$num_states <- sapply(plot_post_tx_log_smoothed_time$term,
function(x){sum(sensitivity_anlys_event_study_data[,x])})
dwplot(plot_post_tx_log_smoothed_time, colour = "black",
vars_order = c(sapply(((max(merged_main_time_data_int$Time_Period_ID) -
min(merged_main_time_data_int$intervention_time_id, na.rm = TRUE)):0),
function(x){paste("pos_", x, "_pd", sep = "")}))) +
theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank(),
panel.background = element_blank(), axis.line = element_line(colour = "black"),
axis.text.x = element_text(angle = 45, size = 4)) +
geom_vline(aes(xintercept = 0), linetype = "dashed") +
labs(y = "Time Periods", x = "Coefficients and 95% Confidence Intervals",
title = "Coefficient of Post-Intervention Periods") +
scale_color_grey() +
coord_flip()
#use this function to compute the cumulative sum, but resets the sum if the variable was 0
compute_cumsum = function(x){
cs = cumsum(x)
cs - cummax((x == 0) * cs)
}
sensitivity_anlys_event_study_data_lin_post_tx <- sensitivity_anlys_event_study_data %>%
arrange(State, Time_Period_ID) %>%
group_by(State) %>%
mutate(
#sum_tx_periods indicates post treatment periods
sum_tx_periods = pos_0_pd + pos_1_pd + pos_2_pd + pos_3_pd +
pos_4_pd + pos_5_pd + pos_6_pd + pos_7_pd + pos_8_pd + pos_9_pd +
pos_10_pd + pos_11_pd + pos_12_pd + pos_13_pd + pos_14_pd +
pos_15_pd + pos_16_pd + pos_17_pd + pos_18_pd + pos_19_pd +
pos_20_pd + pos_21_pd + pos_22_pd + pos_23_pd + pos_24_pd +
pos_25_pd + pos_26_pd + pos_27_pd + pos_28_pd + pos_29_pd +
pos_30_pd + pos_31_pd + pos_32_pd + pos_33_pd + pos_34_pd +
pos_35_pd + pos_36_pd + pos_37_pd + pos_38_pd + pos_39_pd,
#computes the whole number of periods post treatment
time_after_tx = cumsum(sum_tx_periods),
#computes the number of periods post treatment + fractional period when treatment first occurs
num_pd_w_tx = compute_cumsum(Intervention_Redefined ),
#computes the number of periods post the other interventions happening
num_pd_w_naloxone_yes = compute_cumsum(Naloxone_Pharmacy_Yes_Redefined),
num_pd_w_naloxone_no = compute_cumsum(Naloxone_Pharmacy_No_Redefined),
num_pd_w_med_marijuana = compute_cumsum(Medical_Marijuana_Redefined),
num_pd_w_rec_marijuana = compute_cumsum(Recreational_Marijuana_Redefined),
num_pd_w_gsl = compute_cumsum(GSL_Redefined),
num_pd_w_pdmp = compute_cumsum(PDMP_Redefined),
num_pd_w_medicaid = compute_cumsum(Medicaid_Expansion_Redefined),
#compute the lagged versions of the variables so that intercept = effect when tx first occurs
lag_num_pd_w_tx = lag(num_pd_w_tx),
lag_num_pd_w_naloxone_yes = lag(num_pd_w_naloxone_yes),
lag_num_pd_w_naloxone_no = lag(num_pd_w_naloxone_no),
lag_num_pd_w_med_marijuana = lag(num_pd_w_med_marijuana),
lag_num_pd_w_rec_marijuana = lag(num_pd_w_rec_marijuana),
lag_num_pd_w_gsl = lag(num_pd_w_gsl),
lag_num_pd_w_pdmp = lag(num_pd_w_pdmp),
lag_num_pd_w_medicaid = lag(num_pd_w_medicaid))
#fill in a 0 for the NAs so we keep all the data and at most this will be 0
#also fill in a 0 for the linear effect of intervention for the first period
sensitivity_anlys_event_study_data_lin_post_tx$lag_num_pd_w_tx[
sensitivity_anlys_event_study_data_lin_post_tx$Intervention_Redefined < 1] <-
sensitivity_anlys_event_study_data_lin_post_tx$lag_num_pd_w_naloxone_yes[
sensitivity_anlys_event_study_data_lin_post_tx$Naloxone_Pharmacy_Yes_Redefined < 1]<-
sensitivity_anlys_event_study_data_lin_post_tx$lag_num_pd_w_naloxone_no[
sensitivity_anlys_event_study_data_lin_post_tx$Naloxone_Pharmacy_No_Redefined < 1] <-
sensitivity_anlys_event_study_data_lin_post_tx$lag_num_pd_w_med_marijuana[
sensitivity_anlys_event_study_data_lin_post_tx$Medical_Marijuana_Redefined < 1] <-
sensitivity_anlys_event_study_data_lin_post_tx$lag_num_pd_w_rec_marijuana[
sensitivity_anlys_event_study_data_lin_post_tx$Recreational_Marijuana_Redefined < 1] <-
sensitivity_anlys_event_study_data_lin_post_tx$lag_num_pd_w_gsl[
sensitivity_anlys_event_study_data_lin_post_tx$GSL_Redefined < 1] <-
sensitivity_anlys_event_study_data_lin_post_tx$lag_num_pd_w_pdmp[
sensitivity_anlys_event_study_data_lin_post_tx$PDMP_Redefined < 1] <-
sensitivity_anlys_event_study_data_lin_post_tx$lag_num_pd_w_medicaid[
sensitivity_anlys_event_study_data_lin_post_tx$Medicaid_Expansion_Redefined < 1] <-0
#run the gam model
sensitivity_anlys_lin_post_tx_model_log_smoothed_time<-gam(log(prop_dead)~ State +
s(Time_Period_ID, bs = 'cr', by = as.factor(Region)) +
Naloxone_Pharmacy_Yes_Redefined +
lag_num_pd_w_naloxone_yes +
Naloxone_Pharmacy_No_Redefined +
lag_num_pd_w_naloxone_no +
Medical_Marijuana_Redefined +
lag_num_pd_w_med_marijuana +
Recreational_Marijuana_Redefined +
lag_num_pd_w_rec_marijuana +
GSL_Redefined +
lag_num_pd_w_gsl +
PDMP_Redefined +
lag_num_pd_w_pdmp +
Medicaid_Expansion_Redefined +
lag_num_pd_w_medicaid +
Intervention_Redefined +
lag_num_pd_w_tx,
data = sensitivity_anlys_event_study_data_lin_post_tx)
summary(sensitivity_anlys_lin_post_tx_model_log_smoothed_time)
##
## Family: gaussian
## Link function: identity
##
## Formula:
## log(prop_dead) ~ State + s(Time_Period_ID, bs = "cr", by = as.factor(Region)) +
## Naloxone_Pharmacy_Yes_Redefined + lag_num_pd_w_naloxone_yes +
## Naloxone_Pharmacy_No_Redefined + lag_num_pd_w_naloxone_no +
## Medical_Marijuana_Redefined + lag_num_pd_w_med_marijuana +
## Recreational_Marijuana_Redefined + lag_num_pd_w_rec_marijuana +
## GSL_Redefined + lag_num_pd_w_gsl + PDMP_Redefined + lag_num_pd_w_pdmp +
## Medicaid_Expansion_Redefined + lag_num_pd_w_medicaid + Intervention_Redefined +
## lag_num_pd_w_tx
##
## Parametric coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -9.678941 0.062382 -155.157 < 2e-16 ***
## StateAlaska 0.177343 0.083697 2.119 0.034231 *
## StateArizona 0.264886 0.066237 3.999 6.60e-05 ***
## StateArkansas -0.502592 0.067334 -7.464 1.27e-13 ***
## StateCalifornia 0.015999 0.086296 0.185 0.852935
## StateColorado 0.077989 0.078766 0.990 0.322234
## StateConnecticut 0.263003 0.072140 3.646 0.000274 ***
## StateDelaware 0.179210 0.068643 2.611 0.009106 **
## StateFlorida 0.469682 0.068570 6.850 9.97e-12 ***
## StateGeorgia 0.234711 0.070717 3.319 0.000921 ***
## StateHawaii -0.412978 0.086070 -4.798 1.73e-06 ***
## StateIdaho -0.263980 0.070067 -3.768 0.000170 ***
## StateIllinois 0.273254 0.074076 3.689 0.000232 ***
## StateIndiana -0.016469 0.068716 -0.240 0.810617
## StateIowa -0.613142 0.064686 -9.479 < 2e-16 ***
## StateKansas -0.194721 0.064205 -3.033 0.002456 **
## StateKentucky 0.595200 0.069329 8.585 < 2e-16 ***
## StateLouisiana 0.421595 0.064054 6.582 6.00e-11 ***
## StateMaine 0.219451 0.083298 2.635 0.008495 **
## StateMaryland -1.406348 0.070417 -19.972 < 2e-16 ***
## StateMassachusetts -0.143408 0.069758 -2.056 0.039939 *
## StateMichigan 0.008286 0.072491 0.114 0.909009
## StateMinnesota -0.559234 0.070624 -7.918 4.06e-15 ***
## StateMississippi -0.130853 0.063894 -2.048 0.040700 *
## StateMissouri 0.300737 0.073666 4.082 4.64e-05 ***
## StateMontana -0.225004 0.074970 -3.001 0.002724 **
## StateNebraska -0.903276 0.067623 -13.358 < 2e-16 ***
## StateNevada 0.535619 0.081555 6.568 6.59e-11 ***
## StateNew Hampshire 0.177163 0.068298 2.594 0.009561 **
## StateNew Jersey 0.227583 0.067883 3.353 0.000816 ***
## StateNew Mexico 0.684644 0.078563 8.715 < 2e-16 ***
## StateNew York -0.220508 0.070836 -3.113 0.001880 **
## StateNorth Carolina 0.307306 0.063646 4.828 1.49e-06 ***
## StateNorth Dakota -1.208449 0.064240 -18.811 < 2e-16 ***
## StateOhio 0.617865 0.068359 9.039 < 2e-16 ***
## StateOklahoma 0.367308 0.069080 5.317 1.18e-07 ***
## StateOregon -0.092079 0.081816 -1.125 0.260549
## StatePennsylvania 0.639060 0.072385 8.829 < 2e-16 ***
## StateRhode Island -0.312299 0.074771 -4.177 3.09e-05 ***
## StateSouth Carolina 0.126318 0.064639 1.954 0.050825 .
## StateSouth Dakota -1.077450 0.067181 -16.038 < 2e-16 ***
## StateTennessee 0.455862 0.063531 7.175 1.03e-12 ***
## StateTexas 0.065809 0.071916 0.915 0.360265
## StateUtah -0.115434 0.069154 -1.669 0.095236 .
## StateVermont -0.091676 0.070121 -1.307 0.191236
## StateVirginia 0.018175 0.065608 0.277 0.781792
## StateWashington 0.107278 0.080482 1.333 0.182707
## StateWest Virginia 0.656607 0.069869 9.398 < 2e-16 ***
## StateWisconsin 0.106853 0.066196 1.614 0.106655
## StateWyoming -0.069148 0.063912 -1.082 0.279425
## Naloxone_Pharmacy_Yes_Redefined -0.040959 0.041733 -0.981 0.326494
## lag_num_pd_w_naloxone_yes -0.027003 0.006588 -4.099 4.33e-05 ***
## Naloxone_Pharmacy_No_Redefined 0.129651 0.043251 2.998 0.002756 **
## lag_num_pd_w_naloxone_no -0.022777 0.004517 -5.043 5.03e-07 ***
## Medical_Marijuana_Redefined 0.192279 0.030305 6.345 2.78e-10 ***
## lag_num_pd_w_med_marijuana -0.008505 0.002213 -3.843 0.000125 ***
## Recreational_Marijuana_Redefined -0.015978 0.061129 -0.261 0.793829
## lag_num_pd_w_rec_marijuana -0.008431 0.010880 -0.775 0.438514
## GSL_Redefined 0.071508 0.032779 2.182 0.029267 *
## lag_num_pd_w_gsl 0.010557 0.004060 2.600 0.009388 **
## PDMP_Redefined -0.128942 0.027146 -4.750 2.19e-06 ***
## lag_num_pd_w_pdmp 0.006260 0.002417 2.590 0.009666 **
## Medicaid_Expansion_Redefined 0.076951 0.033493 2.298 0.021699 *
## lag_num_pd_w_medicaid 0.012902 0.004989 2.586 0.009776 **
## Intervention_Redefined -0.044881 0.023857 -1.881 0.060090 .
## lag_num_pd_w_tx -0.013751 0.001782 -7.715 1.94e-14 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Approximate significance of smooth terms:
## edf Ref.df F p-value
## s(Time_Period_ID):as.factor(Region)Midwest 4.133 5.113 97.54 <2e-16 ***
## s(Time_Period_ID):as.factor(Region)Northeast 8.398 8.897 49.05 <2e-16 ***
## s(Time_Period_ID):as.factor(Region)South 6.219 7.377 55.25 <2e-16 ***
## s(Time_Period_ID):as.factor(Region)West 4.386 5.424 37.09 <2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## R-sq.(adj) = 0.852 Deviance explained = 85.8%
## GCV = 0.082 Scale est. = 0.078309 n = 1980
plot(sensitivity_anlys_lin_post_tx_model_log_smoothed_time, pages = 1)
#compute the full dataset including basis functions
full_df_w_basis_functions_sensitivity_anlys_lin_post_tx_log_smoothed_time <-
data.frame(predict(sensitivity_anlys_lin_post_tx_model_log_smoothed_time, type = "lpmatrix"))
#estimate the 95% CI and SE
coefficient_values_sensitivity_anlys_lin_post_tx_log_smoothed_time <- coef(sensitivity_anlys_lin_post_tx_model_log_smoothed_time)
sensitivity_anlys_lin_post_tx_sd_and_ci_log_smoothed_time_w_cov <-
compute_sd_and_CI_linear_link(as.matrix(full_df_w_basis_functions_sensitivity_anlys_lin_post_tx_log_smoothed_time),
log(sensitivity_anlys_event_study_data$prop_dead),
coefficient_values_sensitivity_anlys_lin_post_tx_log_smoothed_time,
d = ncol(full_df_w_basis_functions_sensitivity_anlys_lin_post_tx_log_smoothed_time),
return_full_cov = TRUE)
#create dataset for the plot for the coefficients
#since we also returned the variance_covariance matrix here, sensitivity_anlys_lin_post_tx_sd_and_ci_log_smoothed_time_w_cov[[1]]
#contains the table of coefficients
sensitivity_anlys_lin_post_tx_sd_and_ci_log_smoothed_time <- sensitivity_anlys_lin_post_tx_sd_and_ci_log_smoothed_time_w_cov[[1]]
colnames(sensitivity_anlys_lin_post_tx_sd_and_ci_log_smoothed_time) <- c("conf.low", "estimate", "conf.high", "sd")
sensitivity_anlys_lin_post_tx_sd_and_ci_log_smoothed_time$term <- rownames(sensitivity_anlys_lin_post_tx_sd_and_ci_log_smoothed_time)
sensitivity_anlys_lin_post_tx_sd_and_ci_log_smoothed_time$ci_95 <-
paste("95% CI = (", format(round(sensitivity_anlys_lin_post_tx_sd_and_ci_log_smoothed_time$conf.low, 3), nsmall = 3), ", ",
format(round(sensitivity_anlys_lin_post_tx_sd_and_ci_log_smoothed_time$conf.high, 3), nsmall = 3), ")", sep = "")
dwplot(sensitivity_anlys_lin_post_tx_sd_and_ci_log_smoothed_time[51:66,], colour = "black") +
theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank(),
panel.background = element_blank(), axis.line = element_line(colour = "black")) +
geom_vline(aes(xintercept = 0), linetype = "dashed") +
labs(y = "Term", x = "Coefficients and 95% Confidence Intervals",
title = "Coefficient of Analysis With Smoothed Time Effects,
Linear Intervention") +
scale_color_grey() +
geom_text(sensitivity_anlys_lin_post_tx_sd_and_ci_log_smoothed_time[51:66,],
mapping = aes(label = format(round(estimate, 3), nsmall = 3), x = 0.55, y = 16:1), size = 3) +
geom_text(sensitivity_anlys_lin_post_tx_sd_and_ci_log_smoothed_time[51:66,],
mapping = aes(label = ci_95, x = 0.9, y = 16:1), size = 3) +
xlim(-.5, 1.1)
#create table of risk ratios for coefficients
table_of_RR_lin_eff_log_all_yr <- sensitivity_anlys_lin_post_tx_sd_and_ci_log_smoothed_time[51:66,]
table_of_RR_lin_eff_log_all_yr$estimate <- round(exp(table_of_RR_lin_eff_log_all_yr$estimate), 4)
table_of_RR_lin_eff_log_all_yr$conf.low <- exp(table_of_RR_lin_eff_log_all_yr$conf.low)
table_of_RR_lin_eff_log_all_yr$conf.high <- exp(table_of_RR_lin_eff_log_all_yr$conf.high)
table_of_RR_lin_eff_log_all_yr$ci_95 <- paste("95% CI = (",
format(round(table_of_RR_lin_eff_log_all_yr$conf.low, 4), nsmall = 4), ", ",
format(round(table_of_RR_lin_eff_log_all_yr$conf.high, 4), nsmall = 4), ")", sep = "")
# write.csv(table_of_RR_lin_eff_log_all_yr, "./Data/table_of_RR_lin_eff_log_all_yr_5_13_22.csv")
table_of_RR_lin_eff_log_all_yr
## conf.low estimate conf.high sd
## Naloxone_Pharmacy_Yes_Redefined 0.8187863 0.9599 1.1252597 0.081108684
## lag_num_pd_w_naloxone_yes 0.9475055 0.9734 0.9999156 0.013734234
## Naloxone_Pharmacy_No_Redefined 0.9098115 1.1384 1.4244976 0.114371692
## lag_num_pd_w_naloxone_no 0.9590309 0.9775 0.9962848 0.009721913
## Medical_Marijuana_Redefined 1.0112344 1.2120 1.4526453 0.092401666
## lag_num_pd_w_med_marijuana 0.9783527 0.9915 1.0048868 0.006826499
## Recreational_Marijuana_Redefined 0.7995336 0.9841 1.2113932 0.105994347
## lag_num_pd_w_rec_marijuana 0.9538079 0.9916 1.0308994 0.019827701
## GSL_Redefined 0.9396939 1.0741 1.2277916 0.068218936
## lag_num_pd_w_gsl 0.9898149 1.0106 1.0318482 0.010609422
## PDMP_Redefined 0.7515220 0.8790 1.0281596 0.079955398
## lag_num_pd_w_pdmp 0.9893996 1.0063 1.0234486 0.008631352
## Medicaid_Expansion_Redefined 0.9273050 1.0800 1.2578132 0.077767186
## lag_num_pd_w_medicaid 0.9883064 1.0130 1.0382815 0.012584032
## Intervention_Redefined 0.8235221 0.9561 1.1100472 0.076165143
## lag_num_pd_w_tx 0.9734497 0.9863 0.9994067 0.006713178
## term
## Naloxone_Pharmacy_Yes_Redefined Naloxone_Pharmacy_Yes_Redefined
## lag_num_pd_w_naloxone_yes lag_num_pd_w_naloxone_yes
## Naloxone_Pharmacy_No_Redefined Naloxone_Pharmacy_No_Redefined
## lag_num_pd_w_naloxone_no lag_num_pd_w_naloxone_no
## Medical_Marijuana_Redefined Medical_Marijuana_Redefined
## lag_num_pd_w_med_marijuana lag_num_pd_w_med_marijuana
## Recreational_Marijuana_Redefined Recreational_Marijuana_Redefined
## lag_num_pd_w_rec_marijuana lag_num_pd_w_rec_marijuana
## GSL_Redefined GSL_Redefined
## lag_num_pd_w_gsl lag_num_pd_w_gsl
## PDMP_Redefined PDMP_Redefined
## lag_num_pd_w_pdmp lag_num_pd_w_pdmp
## Medicaid_Expansion_Redefined Medicaid_Expansion_Redefined
## lag_num_pd_w_medicaid lag_num_pd_w_medicaid
## Intervention_Redefined Intervention_Redefined
## lag_num_pd_w_tx lag_num_pd_w_tx
## ci_95
## Naloxone_Pharmacy_Yes_Redefined 95% CI = (0.8188, 1.1253)
## lag_num_pd_w_naloxone_yes 95% CI = (0.9475, 0.9999)
## Naloxone_Pharmacy_No_Redefined 95% CI = (0.9098, 1.4245)
## lag_num_pd_w_naloxone_no 95% CI = (0.9590, 0.9963)
## Medical_Marijuana_Redefined 95% CI = (1.0112, 1.4526)
## lag_num_pd_w_med_marijuana 95% CI = (0.9784, 1.0049)
## Recreational_Marijuana_Redefined 95% CI = (0.7995, 1.2114)
## lag_num_pd_w_rec_marijuana 95% CI = (0.9538, 1.0309)
## GSL_Redefined 95% CI = (0.9397, 1.2278)
## lag_num_pd_w_gsl 95% CI = (0.9898, 1.0318)
## PDMP_Redefined 95% CI = (0.7515, 1.0282)
## lag_num_pd_w_pdmp 95% CI = (0.9894, 1.0234)
## Medicaid_Expansion_Redefined 95% CI = (0.9273, 1.2578)
## lag_num_pd_w_medicaid 95% CI = (0.9883, 1.0383)
## Intervention_Redefined 95% CI = (0.8235, 1.1100)
## lag_num_pd_w_tx 95% CI = (0.9734, 0.9994)
date_data <- sensitivity_anlys_event_study_data_lin_post_tx[, c("Time_Period_ID", "Time_Period_Start")]
date_data <- date_data[!duplicated(date_data),]
var_cov_mat_beta_lin_tx <- as.matrix(sensitivity_anlys_lin_post_tx_sd_and_ci_log_smoothed_time_w_cov[[2]][
c("Intervention_Redefined", "lag_num_pd_w_tx"), c("Intervention_Redefined", "lag_num_pd_w_tx")])
attr_deaths_est_log_smoothed_time_lin_post <- attr_death_compute(sensitivity_anlys_event_study_data_lin_post_tx,
sensitivity_anlys_lin_post_tx_sd_and_ci_log_smoothed_time,
tx_name = c("Intervention_Redefined", "lag_num_pd_w_tx"),
var_cov_mat_beta_lin_tx)
attr_deaths_est_log_smoothed_time_lin_post <- merge(attr_deaths_est_log_smoothed_time_lin_post, date_data,
by.x = "Time_Period", by.y = "Time_Period_ID")
attr_deaths_est_log_smoothed_time_lin_post_summary <- attr_deaths_est_log_smoothed_time_lin_post %>%
group_by(year = year(Time_Period_Start)) %>%
summarise(total_attr_deaths = sum(attr_deaths),
total_attr_deaths_lb = sum(attr_deaths_lb),
total_attr_deaths_ub = sum(attr_deaths_ub))
sum(attr_deaths_est_log_smoothed_time_lin_post_summary$total_attr_deaths)
## [1] 197297.8
sum(attr_deaths_est_log_smoothed_time_lin_post_summary$total_attr_deaths_lb)
## [1] -73711.92
sum(attr_deaths_est_log_smoothed_time_lin_post_summary$total_attr_deaths_ub)
## [1] 468307.4
sum(attr_deaths_est_log_smoothed_time_lin_post_summary$total_attr_deaths)/sum(main_analysis_data$imputed_deaths)
## [1] 0.2997463
sum(attr_deaths_est_log_smoothed_time_lin_post_summary$total_attr_deaths_lb)/sum(main_analysis_data$imputed_deaths)
## [1] -0.1119875
sum(attr_deaths_est_log_smoothed_time_lin_post_summary$total_attr_deaths_ub)/sum(main_analysis_data$imputed_deaths)
## [1] 0.71148
ggplot(attr_deaths_est_log_smoothed_time_lin_post_summary, aes(x = year)) +
# geom_point(aes(y = attr_deaths)) +
geom_line(aes(y = total_attr_deaths, linetype = "Estimate")) +
# geom_point(aes(y = attr_deaths_lb)) +
geom_line(aes(y = total_attr_deaths_lb, linetype = "95% CI")) +
# geom_point(aes(y = attr_deaths_ub)) +
geom_line(aes(y = total_attr_deaths_ub, linetype = "95% CI")) +
labs(x = "Date", y = "Yearly Number of Lives Saved",
title = "Estimated Number of Lives Saved Per Year,
Using Smoothed Time Effects,
Log Probability of Drug Overdose Death, Linear Policy Effects",
linetype = "") +
theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank(),
panel.background = element_blank(), axis.line = element_line(colour = "black"),
legend.position = c(.2, .7)) +
scale_linetype_manual(values = c("dashed", "solid"))
#overall national overdose deaths
national_od <- sensitivity_anlys_event_study_data_lin_post_tx %>%
group_by(year = year(Time_Period_Start)) %>%
summarise(total_od = sum(imputed_deaths),
total_od_prob = sum(imputed_deaths)/sum(population),
total_pop = sum(population))
national_od <- merge(national_od, attr_deaths_est_log_smoothed_time_lin_post_summary, by = "year")
ggplot(national_od, aes(x = year)) +
geom_line(aes(y = total_od, color = "Oberved OD")) +
geom_line(aes(y = total_attr_deaths, color = "Lives Saved", linetype = "Estimate")) +
geom_line(aes(y = total_attr_deaths_lb, color = "Lives Saved", linetype = "95% CI")) +
geom_line(aes(y = total_attr_deaths_ub, color = "Lives Saved", linetype = "95% CI")) +
geom_line(aes(y = total_attr_deaths + total_od,
color = "Potential Number of Deaths Had There Not Been DIH Prosecutions")) +
labs(x = "Date", y = "People",
title = "Number of OD Deaths and Estimated Number of Lives Saved in Each Year",
color = "", linetype = "") +
theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank(),
panel.background = element_blank(), axis.line = element_line(colour = "black"),
legend.position = "bottom") +
scale_linetype_manual(values = c("dashed", "solid")) +
guides(color=guide_legend(nrow=2,byrow=TRUE),
linetype=guide_legend(nrow=2,byrow=TRUE))
# pdf("./Figures/lives_saved_ols_lin_tx_all_time_5_13_22.pdf", height = 5, width = 5)
ggplot(attr_deaths_est_log_smoothed_time_lin_post_summary, aes(x = year)) +
# geom_point(aes(y = attr_deaths)) +
geom_line(aes(y = total_attr_deaths, linetype = "Estimate")) +
# geom_point(aes(y = attr_deaths_lb)) +
geom_line(aes(y = total_attr_deaths_lb, linetype = "95% CI")) +
# geom_point(aes(y = attr_deaths_ub)) +
geom_line(aes(y = total_attr_deaths_ub, linetype = "95% CI")) +
labs(x = "Year", y = "Yeary Number of Lives Saved",
# title = "Estimated Number of Lives Saved Per Year,
# Using Smoothed Time Effects,
# Log Probability of Drug Overdose Death, Linear Policy Effects",
linetype = "") +
theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank(),
panel.background = element_blank(), axis.line = element_line(colour = "black"),
legend.position = c(.2, .75)) +
scale_linetype_manual(values = c("dashed", "solid"))
# dev.off()
# pdf("./Figures/lives_saved_and_est_od_ols_lin_tx_all_time_5_13_22.pdf")
ggplot(national_od, aes(x = year)) +
geom_line(aes(y = total_od, color = "Oberved OD")) +
geom_line(aes(y = total_attr_deaths, color = "Lives Saved", linetype = "Estimate")) +
geom_line(aes(y = total_attr_deaths_lb, color = "Lives Saved", linetype = "95% CI")) +
geom_line(aes(y = total_attr_deaths_ub, color = "Lives Saved", linetype = "95% CI")) +
geom_line(aes(y = total_attr_deaths + total_od,
color = "Potential Number of Deaths Had There Not Been DIH Prosecutions")) +
labs(x = "Date", y = "People",
# title = "Number of OD Deaths and Estimated Number of Lives Saved in Each Year",
color = "", linetype = "") +
theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank(),
panel.background = element_blank(), axis.line = element_line(colour = "black"),
legend.position = "bottom") +
scale_linetype_manual(values = c("dashed", "solid")) +
guides(color=guide_legend(nrow=2,byrow=TRUE),
linetype=guide_legend(nrow=2,byrow=TRUE))
# dev.off()
Now we want to run the analysis but on a subset of the periods.
#create a subset dataset where we exclude the last 6 years
data_subset <- sensitivity_anlys_event_study_data_lin_post_tx[sensitivity_anlys_event_study_data_lin_post_tx$Time_Period_ID <= 30,]
sensitivity_anlys_post_tx_model_log_smoothed_time_subset<-gam(log(prop_dead)~ State +
s(Time_Period_ID, bs = 'cr', by = as.factor(Region)) +
Naloxone_Pharmacy_Yes_Redefined +
lag_num_pd_w_naloxone_yes +
Naloxone_Pharmacy_No_Redefined +
lag_num_pd_w_naloxone_no +
Medical_Marijuana_Redefined +
lag_num_pd_w_med_marijuana +
Recreational_Marijuana_Redefined +
lag_num_pd_w_rec_marijuana +
GSL_Redefined +
lag_num_pd_w_gsl +
PDMP_Redefined +
lag_num_pd_w_pdmp +
Medicaid_Expansion_Redefined +
lag_num_pd_w_medicaid +
Intervention_Redefined +
lag_num_pd_w_tx,
data = data_subset)
# summary(sensitivity_anlys_post_tx_model_log_smoothed_time_subset)
#compute the full dataset including basis functions
full_df_w_basis_functions_sensitivity_anlys_post_tx_log_smoothed_time_subset <-
data.frame(predict(sensitivity_anlys_post_tx_model_log_smoothed_time_subset, type = "lpmatrix"))
#estimate the 95% CI and SE
coefficient_values_sensitivity_anlys_post_tx_log_smoothed_time_subset <- coef(sensitivity_anlys_post_tx_model_log_smoothed_time_subset)
sensitivity_anlys_post_tx_sd_and_ci_log_smoothed_time_subset_w_cov <-
compute_sd_and_CI_linear_link(as.matrix(full_df_w_basis_functions_sensitivity_anlys_post_tx_log_smoothed_time_subset),
log(data_subset$prop_dead),
coefficient_values_sensitivity_anlys_post_tx_log_smoothed_time_subset,
d = ncol(full_df_w_basis_functions_sensitivity_anlys_post_tx_log_smoothed_time_subset),
return_full_cov = TRUE)
# sensitivity_anlys_post_tx_sd_and_ci_log_smoothed_time_subset
#since we also returned the variance_covariance matrix here, sensitivity_anlys_post_tx_sd_and_ci_log_smoothed_time_subset_w_cov[[1]]
sensitivity_anlys_post_tx_sd_and_ci_log_smoothed_time_subset <- sensitivity_anlys_post_tx_sd_and_ci_log_smoothed_time_subset_w_cov[[1]]
colnames(sensitivity_anlys_post_tx_sd_and_ci_log_smoothed_time_subset) <- c("conf.low", "estimate", "conf.high", "sd")
sensitivity_anlys_post_tx_sd_and_ci_log_smoothed_time_subset$term <- rownames(sensitivity_anlys_post_tx_sd_and_ci_log_smoothed_time_subset)
sensitivity_anlys_post_tx_sd_and_ci_log_smoothed_time_subset$ci_95 <-
paste("95% CI = (", format(round(sensitivity_anlys_post_tx_sd_and_ci_log_smoothed_time_subset$conf.low, 3), nsmall = 3), ", ",
format(round(sensitivity_anlys_post_tx_sd_and_ci_log_smoothed_time_subset$conf.high, 3), nsmall = 3), ")", sep = "")
dwplot(sensitivity_anlys_post_tx_sd_and_ci_log_smoothed_time_subset[51:66,], colour = "black") +
theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank(),
panel.background = element_blank(), axis.line = element_line(colour = "black")) +
geom_vline(aes(xintercept = 0), linetype = "dashed") +
labs(y = "Term", x = "Coefficients and 95% Confidence Intervals",
title = "Coefficient of Analysis With Smoothed Time Effects,
Linear Intervention, Subset Data") +
scale_color_grey() +
geom_text(sensitivity_anlys_post_tx_sd_and_ci_log_smoothed_time_subset[51:66,],
mapping = aes(label = format(round(estimate, 3), nsmall = 3), x = 0.55, y = 16:1), size = 3) +
geom_text(sensitivity_anlys_post_tx_sd_and_ci_log_smoothed_time_subset[51:66,],
mapping = aes(label = ci_95, x = 0.9, y = 16:1), size = 3) +
xlim(-.7, 1.1)
table_of_RR_lin_eff_log_all_yr_subset <- sensitivity_anlys_post_tx_sd_and_ci_log_smoothed_time_subset[51:66,]
table_of_RR_lin_eff_log_all_yr_subset$estimate <- round(exp(table_of_RR_lin_eff_log_all_yr_subset$estimate), 3)
table_of_RR_lin_eff_log_all_yr_subset$conf.low <- exp(table_of_RR_lin_eff_log_all_yr_subset$conf.low)
table_of_RR_lin_eff_log_all_yr_subset$conf.high <- exp(table_of_RR_lin_eff_log_all_yr_subset$conf.high)
table_of_RR_lin_eff_log_all_yr_subset$ci_95 <- paste("95% CI = (",
format(round(table_of_RR_lin_eff_log_all_yr_subset$conf.low, 3), nsmall = 3), ", ",
format(round(table_of_RR_lin_eff_log_all_yr_subset$conf.high, 3), nsmall = 3), ")",
sep = "")
# write.csv(table_of_RR_lin_eff_log_all_yr_subset, "./Data/table_of_RR_lin_eff_log_all_yr_subset_5_13_22.csv")
table_of_RR_lin_eff_log_all_yr_subset
## conf.low estimate conf.high sd
## Naloxone_Pharmacy_Yes_Redefined 0.6631265 0.836 1.0551939 0.118498452
## lag_num_pd_w_naloxone_yes 0.8849338 0.943 1.0050832 0.032477728
## Naloxone_Pharmacy_No_Redefined 0.8467222 1.162 1.5951617 0.161570850
## lag_num_pd_w_naloxone_no 0.9322543 0.957 0.9818023 0.013210280
## Medical_Marijuana_Redefined 0.9875543 1.264 1.6187505 0.126065909
## lag_num_pd_w_med_marijuana 0.9678008 0.984 1.0007205 0.008532959
## Recreational_Marijuana_Redefined 0.5177888 0.779 1.1723436 0.208467505
## lag_num_pd_w_rec_marijuana 0.8761929 1.015 1.1753395 0.074930112
## GSL_Redefined 0.9001522 1.144 1.4545583 0.122421870
## lag_num_pd_w_gsl 0.9714048 1.010 1.0495745 0.019744107
## PDMP_Redefined 0.7548766 0.884 1.0358000 0.080707939
## lag_num_pd_w_pdmp 0.9925185 1.011 1.0298725 0.009424631
## Medicaid_Expansion_Redefined 0.8467866 1.034 1.2637889 0.102148186
## lag_num_pd_w_medicaid 0.9574030 1.016 1.0786728 0.030424048
## Intervention_Redefined 0.8428284 0.995 1.1746876 0.084692380
## lag_num_pd_w_tx 0.9694544 0.983 0.9971206 0.007178135
## term
## Naloxone_Pharmacy_Yes_Redefined Naloxone_Pharmacy_Yes_Redefined
## lag_num_pd_w_naloxone_yes lag_num_pd_w_naloxone_yes
## Naloxone_Pharmacy_No_Redefined Naloxone_Pharmacy_No_Redefined
## lag_num_pd_w_naloxone_no lag_num_pd_w_naloxone_no
## Medical_Marijuana_Redefined Medical_Marijuana_Redefined
## lag_num_pd_w_med_marijuana lag_num_pd_w_med_marijuana
## Recreational_Marijuana_Redefined Recreational_Marijuana_Redefined
## lag_num_pd_w_rec_marijuana lag_num_pd_w_rec_marijuana
## GSL_Redefined GSL_Redefined
## lag_num_pd_w_gsl lag_num_pd_w_gsl
## PDMP_Redefined PDMP_Redefined
## lag_num_pd_w_pdmp lag_num_pd_w_pdmp
## Medicaid_Expansion_Redefined Medicaid_Expansion_Redefined
## lag_num_pd_w_medicaid lag_num_pd_w_medicaid
## Intervention_Redefined Intervention_Redefined
## lag_num_pd_w_tx lag_num_pd_w_tx
## ci_95
## Naloxone_Pharmacy_Yes_Redefined 95% CI = (0.663, 1.055)
## lag_num_pd_w_naloxone_yes 95% CI = (0.885, 1.005)
## Naloxone_Pharmacy_No_Redefined 95% CI = (0.847, 1.595)
## lag_num_pd_w_naloxone_no 95% CI = (0.932, 0.982)
## Medical_Marijuana_Redefined 95% CI = (0.988, 1.619)
## lag_num_pd_w_med_marijuana 95% CI = (0.968, 1.001)
## Recreational_Marijuana_Redefined 95% CI = (0.518, 1.172)
## lag_num_pd_w_rec_marijuana 95% CI = (0.876, 1.175)
## GSL_Redefined 95% CI = (0.900, 1.455)
## lag_num_pd_w_gsl 95% CI = (0.971, 1.050)
## PDMP_Redefined 95% CI = (0.755, 1.036)
## lag_num_pd_w_pdmp 95% CI = (0.993, 1.030)
## Medicaid_Expansion_Redefined 95% CI = (0.847, 1.264)
## lag_num_pd_w_medicaid 95% CI = (0.957, 1.079)
## Intervention_Redefined 95% CI = (0.843, 1.175)
## lag_num_pd_w_tx 95% CI = (0.969, 0.997)
date_data_subset <- data_subset[, c("Time_Period_ID", "Time_Period_Start")]
date_data_subset <- date_data_subset[!duplicated(date_data_subset),]
var_cov_mat_beta_lin_tx_subset <- as.matrix(sensitivity_anlys_post_tx_sd_and_ci_log_smoothed_time_subset_w_cov[[2]][
c("Intervention_Redefined", "lag_num_pd_w_tx"), c("Intervention_Redefined", "lag_num_pd_w_tx")])
attr_deaths_est_log_smoothed_time_lin_post_subset <- attr_death_compute(data_subset,
sensitivity_anlys_post_tx_sd_and_ci_log_smoothed_time_subset,
tx_name = c("Intervention_Redefined", "lag_num_pd_w_tx"),
var_cov_mat_beta_lin_tx_subset)
attr_deaths_est_log_smoothed_time_lin_post_subset <- merge(attr_deaths_est_log_smoothed_time_lin_post_subset, date_data_subset,
by.x = "Time_Period", by.y = "Time_Period_ID")
attr_deaths_est_log_smoothed_time_lin_post_subset_summary <- attr_deaths_est_log_smoothed_time_lin_post_subset %>%
group_by(year = year(Time_Period_Start)) %>%
summarise(total_attr_deaths = sum(attr_deaths),
total_attr_deaths_lb = sum(attr_deaths_lb),
total_attr_deaths_ub = sum(attr_deaths_ub))
sum(attr_deaths_est_log_smoothed_time_lin_post_subset_summary$total_attr_deaths)
## [1] 67149.1
sum(attr_deaths_est_log_smoothed_time_lin_post_subset_summary$total_attr_deaths_lb)
## [1] -31463.76
sum(attr_deaths_est_log_smoothed_time_lin_post_subset_summary$total_attr_deaths_ub)
## [1] 165762
sum(attr_deaths_est_log_smoothed_time_lin_post_subset_summary$total_attr_deaths)/sum(data_subset$imputed_deaths)
## [1] 0.1771394
sum(attr_deaths_est_log_smoothed_time_lin_post_subset_summary$total_attr_deaths_lb)/sum(data_subset$imputed_deaths)
## [1] -0.08300144
sum(attr_deaths_est_log_smoothed_time_lin_post_subset_summary$total_attr_deaths_ub)/sum(data_subset$imputed_deaths)
## [1] 0.4372803
ggplot(attr_deaths_est_log_smoothed_time_lin_post_subset_summary, aes(x = year)) +
# geom_point(aes(y = attr_deaths)) +
geom_line(aes(y = total_attr_deaths, linetype = "Estimate")) +
# geom_point(aes(y = attr_deaths_lb)) +
geom_line(aes(y = total_attr_deaths_lb, linetype = "95% CI")) +
# geom_point(aes(y = attr_deaths_ub)) +
geom_line(aes(y = total_attr_deaths_ub, linetype = "95% CI")) +
labs(x = "Date", y = "Lives Saved",
title = "Estimated Number of Lives Saved Per Year
Using Smoothed Time Effects,
Log Probability of Drug Overdose Death, Linear Policy Effects",
linetype = "") +
theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank(),
panel.background = element_blank(), axis.line = element_line(colour = "black")) +
scale_linetype_manual(values = c("dashed", "solid"))
ggplot() +
geom_line(attr_deaths_est_log_smoothed_time_lin_post_summary,
mapping = aes(x = year, y = total_attr_deaths, linetype = "Estimate",
color = "Full Dataset")) +
geom_line(attr_deaths_est_log_smoothed_time_lin_post_summary,
mapping = aes(x = year, y = total_attr_deaths_lb, linetype = "95% CI",
color = "Full Dataset")) +
geom_line(attr_deaths_est_log_smoothed_time_lin_post_summary,
mapping = aes(x = year, y = total_attr_deaths_ub, linetype = "95% CI",
color = "Full Dataset")) +
geom_line(attr_deaths_est_log_smoothed_time_lin_post_subset_summary,
mapping = aes(x = year, y = total_attr_deaths, linetype = "Estimate",
color = "Subset Dataset")) +
geom_line(attr_deaths_est_log_smoothed_time_lin_post_subset_summary,
mapping = aes(x = year, y = total_attr_deaths_lb,
linetype = "95% CI",
color = "Subset Dataset")) +
geom_line(attr_deaths_est_log_smoothed_time_lin_post_subset_summary,
mapping = aes(x = year, y = total_attr_deaths_ub,
linetype = "95% CI",
color = "Subset Dataset")) +
labs(x = "Date", y = "Lives Saved",
title = "Estimated Number of Lives Saved Using Smoothed Time Effects,
Log Probability of Drug Overdose Death, Linear Policy Effects",
linetype = "", color = "") +
theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank(),
panel.background = element_blank(), axis.line = element_line(colour = "black")) +
scale_linetype_manual(values = c("dashed", "solid"))
#run the gam model
sensitivity_anlys_lin_post_tx_model_logistic_smoothed_time<-gam(cbind(round(imputed_deaths), round(num_alive))~ State +
s(Time_Period_ID, bs = 'cr', by = as.factor(Region)) +
Naloxone_Pharmacy_Yes_Redefined +
lag_num_pd_w_naloxone_yes +
Naloxone_Pharmacy_No_Redefined +
lag_num_pd_w_naloxone_no +
Medical_Marijuana_Redefined +
lag_num_pd_w_med_marijuana +
Recreational_Marijuana_Redefined +
lag_num_pd_w_rec_marijuana +
GSL_Redefined +
lag_num_pd_w_gsl +
PDMP_Redefined +
lag_num_pd_w_pdmp +
Medicaid_Expansion_Redefined +
lag_num_pd_w_medicaid +
Intervention_Redefined +
lag_num_pd_w_tx,
data = sensitivity_anlys_event_study_data_lin_post_tx,
family = "binomial")
summary(sensitivity_anlys_lin_post_tx_model_logistic_smoothed_time)
##
## Family: binomial
## Link function: logit
##
## Formula:
## cbind(round(imputed_deaths), round(num_alive)) ~ State + s(Time_Period_ID,
## bs = "cr", by = as.factor(Region)) + Naloxone_Pharmacy_Yes_Redefined +
## lag_num_pd_w_naloxone_yes + Naloxone_Pharmacy_No_Redefined +
## lag_num_pd_w_naloxone_no + Medical_Marijuana_Redefined +
## lag_num_pd_w_med_marijuana + Recreational_Marijuana_Redefined +
## lag_num_pd_w_rec_marijuana + GSL_Redefined + lag_num_pd_w_gsl +
## PDMP_Redefined + lag_num_pd_w_pdmp + Medicaid_Expansion_Redefined +
## lag_num_pd_w_medicaid + Intervention_Redefined + lag_num_pd_w_tx
##
## Parametric coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -9.6918398 0.0165286 -586.368 < 2e-16 ***
## StateAlaska 0.1660062 0.0316907 5.238 1.62e-07 ***
## StateArizona 0.2470040 0.0148999 16.578 < 2e-16 ***
## StateArkansas -0.4514581 0.0211955 -21.300 < 2e-16 ***
## StateCalifornia 0.0115038 0.0191547 0.601 0.548126
## StateColorado 0.1212188 0.0196259 6.176 6.56e-10 ***
## StateConnecticut 0.0859484 0.0173084 4.966 6.85e-07 ***
## StateDelaware 0.2921966 0.0245539 11.900 < 2e-16 ***
## StateFlorida 0.4469502 0.0148057 30.188 < 2e-16 ***
## StateGeorgia 0.2005269 0.0169730 11.814 < 2e-16 ***
## StateHawaii -0.3867800 0.0302216 -12.798 < 2e-16 ***
## StateIdaho -0.3168152 0.0254657 -12.441 < 2e-16 ***
## StateIllinois 0.1658041 0.0165498 10.018 < 2e-16 ***
## StateIndiana -0.0491962 0.0157645 -3.121 0.001804 **
## StateIowa -0.6671034 0.0219000 -30.461 < 2e-16 ***
## StateKansas -0.3062736 0.0197456 -15.511 < 2e-16 ***
## StateKentucky 0.5188091 0.0154340 33.615 < 2e-16 ***
## StateLouisiana 0.4010566 0.0149643 26.801 < 2e-16 ***
## StateMaine 0.2586864 0.0263564 9.815 < 2e-16 ***
## StateMaryland -0.9568958 0.0217895 -43.915 < 2e-16 ***
## StateMassachusetts 0.1586346 0.0151818 10.449 < 2e-16 ***
## StateMichigan 0.0130877 0.0162585 0.805 0.420834
## StateMinnesota -0.5908795 0.0197106 -29.978 < 2e-16 ***
## StateMississippi -0.2354073 0.0188147 -12.512 < 2e-16 ***
## StateMissouri 0.3211547 0.0188890 17.002 < 2e-16 ***
## StateMontana -0.1494454 0.0313434 -4.768 1.86e-06 ***
## StateNebraska -0.9918381 0.0296070 -33.500 < 2e-16 ***
## StateNevada 0.5342117 0.0205026 26.056 < 2e-16 ***
## StateNew Hampshire 0.2682052 0.0215241 12.461 < 2e-16 ***
## StateNew Jersey 0.2753953 0.0154320 17.846 < 2e-16 ***
## StateNew Mexico 0.5017865 0.0195363 25.685 < 2e-16 ***
## StateNew York -0.3234978 0.0143800 -22.496 < 2e-16 ***
## StateNorth Carolina 0.2718936 0.0132117 20.580 < 2e-16 ***
## StateNorth Dakota -1.1953770 0.0457733 -26.115 < 2e-16 ***
## StateOhio 0.6679880 0.0152538 43.792 < 2e-16 ***
## StateOklahoma 0.2955589 0.0165172 17.894 < 2e-16 ***
## StateOregon -0.0682414 0.0223395 -3.055 0.002253 **
## StatePennsylvania 0.5992187 0.0159711 37.519 < 2e-16 ***
## StateRhode Island 0.1401910 0.0245331 5.714 1.10e-08 ***
## StateSouth Carolina 0.0672272 0.0160186 4.197 2.71e-05 ***
## StateSouth Dakota -1.1228077 0.0435739 -25.768 < 2e-16 ***
## StateTennessee 0.4346807 0.0133160 32.643 < 2e-16 ***
## StateTexas -0.0454967 0.0159415 -2.854 0.004317 **
## StateUtah 0.0780823 0.0196717 3.969 7.21e-05 ***
## StateVermont -0.1154520 0.0325629 -3.546 0.000392 ***
## StateVirginia -0.0106755 0.0150892 -0.707 0.479259
## StateWashington 0.1034138 0.0192418 5.374 7.68e-08 ***
## StateWest Virginia 0.6933803 0.0173454 39.975 < 2e-16 ***
## StateWisconsin 0.0702017 0.0158598 4.426 9.58e-06 ***
## StateWyoming -0.0314449 0.0337961 -0.930 0.352148
## Naloxone_Pharmacy_Yes_Redefined 0.0155747 0.0080655 1.931 0.053480 .
## lag_num_pd_w_naloxone_yes -0.0213838 0.0011158 -19.165 < 2e-16 ***
## Naloxone_Pharmacy_No_Redefined 0.0619851 0.0079301 7.816 5.43e-15 ***
## lag_num_pd_w_naloxone_no -0.0021388 0.0009269 -2.307 0.021031 *
## Medical_Marijuana_Redefined 0.0335407 0.0058322 5.751 8.87e-09 ***
## lag_num_pd_w_med_marijuana -0.0038222 0.0006421 -5.953 2.64e-09 ***
## Recreational_Marijuana_Redefined 0.0018514 0.0101063 0.183 0.854649
## lag_num_pd_w_rec_marijuana -0.0199844 0.0019355 -10.325 < 2e-16 ***
## GSL_Redefined 0.0256909 0.0060998 4.212 2.53e-05 ***
## lag_num_pd_w_gsl 0.0135896 0.0008501 15.987 < 2e-16 ***
## PDMP_Redefined -0.0327741 0.0068817 -4.762 1.91e-06 ***
## lag_num_pd_w_pdmp 0.0046873 0.0006373 7.355 1.91e-13 ***
## Medicaid_Expansion_Redefined 0.0676747 0.0067899 9.967 < 2e-16 ***
## lag_num_pd_w_medicaid 0.0142307 0.0010244 13.892 < 2e-16 ***
## Intervention_Redefined -0.0187494 0.0063105 -2.971 0.002967 **
## lag_num_pd_w_tx -0.0148682 0.0005601 -26.544 < 2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Approximate significance of smooth terms:
## edf Ref.df Chi.sq p-value
## s(Time_Period_ID):as.factor(Region)Midwest 8.923 8.998 4729 <2e-16 ***
## s(Time_Period_ID):as.factor(Region)Northeast 8.964 9.000 4285 <2e-16 ***
## s(Time_Period_ID):as.factor(Region)South 8.937 8.999 4379 <2e-16 ***
## s(Time_Period_ID):as.factor(Region)West 8.775 8.984 2117 <2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## R-sq.(adj) = 0.917 Deviance explained = 91.2%
## UBRE = 7.8233 Scale est. = 1 n = 1980
plot(sensitivity_anlys_lin_post_tx_model_logistic_smoothed_time, pages = 1)
t#compute the full dataset including basis functions
## standardGeneric for "t" defined from package "base"
##
## function (x)
## standardGeneric("t")
## <environment: 0x7fe0d546dee8>
## Methods may be defined for arguments: x
## Use showMethods("t") for currently available ones.
full_df_w_basis_functions_sensitivity_anlys_lin_post_tx_logistic_smoothed_time <-
data.frame(predict(sensitivity_anlys_lin_post_tx_model_logistic_smoothed_time, type = "lpmatrix"))
#estimate the 95% CI and SD
coefficient_values_sensitivity_anlys_lin_post_tx_logistic_smoothed_time <- coef(sensitivity_anlys_lin_post_tx_model_logistic_smoothed_time)
sensitivity_anlys_lin_post_tx_sd_and_ci_logistic_smoothed_time_w_cov <-
compute_sd_and_CI_logistic_link(as.matrix(full_df_w_basis_functions_sensitivity_anlys_lin_post_tx_logistic_smoothed_time),
(sensitivity_anlys_event_study_data$prop_dead),
coefficient_values_sensitivity_anlys_lin_post_tx_logistic_smoothed_time,
d = ncol(full_df_w_basis_functions_sensitivity_anlys_lin_post_tx_logistic_smoothed_time),
return_full_cov = TRUE)
#create the dataset for the plot
sensitivity_anlys_lin_post_tx_sd_and_ci_logistic_smoothed_time <- sensitivity_anlys_lin_post_tx_sd_and_ci_logistic_smoothed_time_w_cov[[1]]
colnames(sensitivity_anlys_lin_post_tx_sd_and_ci_logistic_smoothed_time) <- c("conf.low", "estimate", "conf.high", "exp_lb",
"exp_coef", "exp_ub", "sd")
sensitivity_anlys_lin_post_tx_sd_and_ci_logistic_smoothed_time$term <- rownames(sensitivity_anlys_lin_post_tx_sd_and_ci_logistic_smoothed_time)
sensitivity_anlys_lin_post_tx_sd_and_ci_logistic_smoothed_time$ci_95 <-
paste("95% CI = (", format(round(sensitivity_anlys_lin_post_tx_sd_and_ci_logistic_smoothed_time$conf.low, 3), nsmall = 3), ", ",
format(round(sensitivity_anlys_lin_post_tx_sd_and_ci_logistic_smoothed_time$conf.high, 3), nsmall = 3), ")", sep = "")
dwplot(sensitivity_anlys_lin_post_tx_sd_and_ci_logistic_smoothed_time[51:66,], colour = "black") +
theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank(),
panel.background = element_blank(), axis.line = element_line(colour = "black")) +
geom_vline(aes(xintercept = 0), linetype = "dashed") +
labs(y = "Term", x = "Coefficients and 95% Confidence Intervals",
title = "Coefficient of Analysis With Smoothed Time Effects,
Linear Intervention") +
scale_color_grey() +
geom_text(sensitivity_anlys_lin_post_tx_sd_and_ci_logistic_smoothed_time[51:66,],
mapping = aes(label = format(round(estimate, 3), nsmall = 3), x = 0.55, y = 16:1), size = 3) +
geom_text(sensitivity_anlys_lin_post_tx_sd_and_ci_logistic_smoothed_time[51:66,],
mapping = aes(label = ci_95, x = 0.9, y = 16:1), size = 3) +
xlim(-.5, 1.1)
table_of_RR_lin_eff_logistic_all_yr <- sensitivity_anlys_lin_post_tx_sd_and_ci_logistic_smoothed_time[51:66,]
table_of_RR_lin_eff_logistic_all_yr$estimate <- round(exp(table_of_RR_lin_eff_logistic_all_yr$estimate), 3)
table_of_RR_lin_eff_logistic_all_yr$conf.low <- exp(table_of_RR_lin_eff_logistic_all_yr$conf.low)
table_of_RR_lin_eff_logistic_all_yr$conf.high <- exp(table_of_RR_lin_eff_logistic_all_yr$conf.high)
table_of_RR_lin_eff_logistic_all_yr$ci_95 <- paste("95% CI = (",
format(round(table_of_RR_lin_eff_logistic_all_yr$conf.low, 3), nsmall = 3), ", ",
format(round(table_of_RR_lin_eff_logistic_all_yr$conf.high, 3), nsmall = 3), ")",
sep = "")
# write.csv(table_of_RR_lin_eff_logistic_all_yr, "./Data/table_of_RR_lin_eff_logistic_all_yr_5_13_22.csv")
table_of_RR_lin_eff_logistic_all_yr
## conf.low estimate conf.high exp_lb
## Naloxone_Pharmacy_Yes_Redefined 0.8731089 1.016 1.181570 0.8731089
## lag_num_pd_w_naloxone_yes 0.9561626 0.979 1.002062 0.9561626
## Naloxone_Pharmacy_No_Redefined 0.8912458 1.064 1.270112 0.8912458
## lag_num_pd_w_naloxone_no 0.9810328 0.998 1.014983 0.9810328
## Medical_Marijuana_Redefined 0.8811402 1.034 1.213635 0.8811402
## lag_num_pd_w_med_marijuana 0.9806591 0.996 1.011957 0.9806591
## Recreational_Marijuana_Redefined 0.8578127 1.002 1.170080 0.8578127
## lag_num_pd_w_rec_marijuana 0.9493872 0.980 1.012042 0.9493872
## GSL_Redefined 0.9106076 1.026 1.156069 0.9106076
## lag_num_pd_w_gsl 0.9935669 1.014 1.034205 0.9935669
## PDMP_Redefined 0.8171409 0.968 1.146135 0.8171409
## lag_num_pd_w_pdmp 0.9832284 1.005 1.026637 0.9832284
## Medicaid_Expansion_Redefined 0.9529083 1.070 1.201518 0.9529083
## lag_num_pd_w_medicaid 0.9949598 1.014 1.034082 0.9949598
## Intervention_Redefined 0.8347624 0.981 1.153856 0.8347624
## lag_num_pd_w_tx 0.9667396 0.985 1.004098 0.9667396
## exp_coef exp_ub sd
## Naloxone_Pharmacy_Yes_Redefined 1.0156966 1.181570 0.077178405
## lag_num_pd_w_naloxone_yes 0.9788432 1.002062 0.011960947
## Naloxone_Pharmacy_No_Redefined 1.0639465 1.270112 0.090367361
## lag_num_pd_w_naloxone_no 0.9978635 1.014983 0.008678879
## Medical_Marijuana_Redefined 1.0341096 1.213635 0.081673107
## lag_num_pd_w_med_marijuana 0.9961851 1.011957 0.008014361
## Recreational_Marijuana_Redefined 1.0018531 1.170080 0.079194291
## lag_num_pd_w_rec_marijuana 0.9802140 1.012042 0.016303184
## GSL_Redefined 1.0260238 1.156069 0.060884760
## lag_num_pd_w_gsl 1.0136824 1.034205 0.010226281
## PDMP_Redefined 0.9677572 1.146135 0.086311045
## lag_num_pd_w_pdmp 1.0046983 1.026637 0.011021019
## Medicaid_Expansion_Redefined 1.0700172 1.201518 0.059138447
## lag_num_pd_w_medicaid 1.0143325 1.034082 0.009838603
## Intervention_Redefined 0.9814253 1.153856 0.082581000
## lag_num_pd_w_tx 0.9852418 1.004098 0.009672384
## term
## Naloxone_Pharmacy_Yes_Redefined Naloxone_Pharmacy_Yes_Redefined
## lag_num_pd_w_naloxone_yes lag_num_pd_w_naloxone_yes
## Naloxone_Pharmacy_No_Redefined Naloxone_Pharmacy_No_Redefined
## lag_num_pd_w_naloxone_no lag_num_pd_w_naloxone_no
## Medical_Marijuana_Redefined Medical_Marijuana_Redefined
## lag_num_pd_w_med_marijuana lag_num_pd_w_med_marijuana
## Recreational_Marijuana_Redefined Recreational_Marijuana_Redefined
## lag_num_pd_w_rec_marijuana lag_num_pd_w_rec_marijuana
## GSL_Redefined GSL_Redefined
## lag_num_pd_w_gsl lag_num_pd_w_gsl
## PDMP_Redefined PDMP_Redefined
## lag_num_pd_w_pdmp lag_num_pd_w_pdmp
## Medicaid_Expansion_Redefined Medicaid_Expansion_Redefined
## lag_num_pd_w_medicaid lag_num_pd_w_medicaid
## Intervention_Redefined Intervention_Redefined
## lag_num_pd_w_tx lag_num_pd_w_tx
## ci_95
## Naloxone_Pharmacy_Yes_Redefined 95% CI = (0.873, 1.182)
## lag_num_pd_w_naloxone_yes 95% CI = (0.956, 1.002)
## Naloxone_Pharmacy_No_Redefined 95% CI = (0.891, 1.270)
## lag_num_pd_w_naloxone_no 95% CI = (0.981, 1.015)
## Medical_Marijuana_Redefined 95% CI = (0.881, 1.214)
## lag_num_pd_w_med_marijuana 95% CI = (0.981, 1.012)
## Recreational_Marijuana_Redefined 95% CI = (0.858, 1.170)
## lag_num_pd_w_rec_marijuana 95% CI = (0.949, 1.012)
## GSL_Redefined 95% CI = (0.911, 1.156)
## lag_num_pd_w_gsl 95% CI = (0.994, 1.034)
## PDMP_Redefined 95% CI = (0.817, 1.146)
## lag_num_pd_w_pdmp 95% CI = (0.983, 1.027)
## Medicaid_Expansion_Redefined 95% CI = (0.953, 1.202)
## lag_num_pd_w_medicaid 95% CI = (0.995, 1.034)
## Intervention_Redefined 95% CI = (0.835, 1.154)
## lag_num_pd_w_tx 95% CI = (0.967, 1.004)
date_data <- sensitivity_anlys_event_study_data_lin_post_tx[, c("Time_Period_ID", "Time_Period_Start")]
date_data <- date_data[!duplicated(date_data),]
var_cov_mat_beta_logistic_tx <- as.matrix(sensitivity_anlys_lin_post_tx_sd_and_ci_logistic_smoothed_time_w_cov[[2]][
c("Intervention_Redefined", "lag_num_pd_w_tx"), c("Intervention_Redefined", "lag_num_pd_w_tx")])
attr_deaths_est_logistic_smoothed_time_lin_post <- attr_death_compute(sensitivity_anlys_event_study_data_lin_post_tx,
sensitivity_anlys_lin_post_tx_sd_and_ci_logistic_smoothed_time,
tx_name = c("Intervention_Redefined",
"lag_num_pd_w_tx"),
var_cov_mat_beta_logistic_tx)
attr_deaths_est_logistic_smoothed_time_lin_post <- merge(attr_deaths_est_logistic_smoothed_time_lin_post, date_data,
by.x = "Time_Period", by.y = "Time_Period_ID")
attr_deaths_est_logistic_smoothed_time_lin_post_summary <- attr_deaths_est_logistic_smoothed_time_lin_post %>%
group_by(year = year(Time_Period_Start)) %>%
summarise(total_attr_deaths = sum(attr_deaths),
total_attr_deaths_lb = sum(attr_deaths_lb),
total_attr_deaths_ub = sum(attr_deaths_ub))
sum(attr_deaths_est_logistic_smoothed_time_lin_post_summary$total_attr_deaths)
## [1] 194246
sum(attr_deaths_est_logistic_smoothed_time_lin_post_summary$total_attr_deaths_lb)
## [1] -179850.4
sum(attr_deaths_est_logistic_smoothed_time_lin_post_summary$total_attr_deaths_ub)
## [1] 568342.3
sum(attr_deaths_est_logistic_smoothed_time_lin_post_summary$total_attr_deaths)/sum(main_analysis_data$imputed_deaths)
## [1] 0.2951098
sum(attr_deaths_est_logistic_smoothed_time_lin_post_summary$total_attr_deaths_lb)/sum(main_analysis_data$imputed_deaths)
## [1] -0.2732392
sum(attr_deaths_est_logistic_smoothed_time_lin_post_summary$total_attr_deaths_ub)/sum(main_analysis_data$imputed_deaths)
## [1] 0.8634589
ggplot(attr_deaths_est_logistic_smoothed_time_lin_post_summary, aes(x = year)) +
# geom_point(aes(y = attr_deaths)) +
geom_line(aes(y = total_attr_deaths, linetype = "Estimate")) +
# geom_point(aes(y = attr_deaths_lb)) +
geom_line(aes(y = total_attr_deaths_lb, linetype = "95% CI")) +
# geom_point(aes(y = attr_deaths_ub)) +
geom_line(aes(y = total_attr_deaths_ub, linetype = "95% CI")) +
labs(x = "Date", y = "Lives Saved",
title = "Estimated Number of Lives Saved Per Year,
Using Smoothed Time Effects,
Log Probability of Drug Overdose Death, Linear Policy Effects",
linetype = "") +
theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank(),
panel.background = element_blank(), axis.line = element_line(colour = "black")) +
scale_linetype_manual(values = c("dashed", "solid"))
#overall national overdose deaths
national_od <- sensitivity_anlys_event_study_data_lin_post_tx %>%
group_by(year = year(Time_Period_Start)) %>%
summarise(total_od = sum(imputed_deaths),
total_od_prob = sum(imputed_deaths)/sum(population),
total_pop = sum(population))
national_od <- merge(national_od, attr_deaths_est_logistic_smoothed_time_lin_post_summary, by = "year")
ggplot(national_od, aes(x = year)) +
geom_line(aes(y = total_od, color = "Oberved OD")) +
geom_line(aes(y = total_attr_deaths, color = "Lives Saved", linetype = "Estimate")) +
geom_line(aes(y = total_attr_deaths_lb, color = "Lives Saved", linetype = "95% CI")) +
geom_line(aes(y = total_attr_deaths_ub, color = "Lives Saved", linetype = "95% CI")) +
geom_line(aes(y = total_attr_deaths + total_od,
color = "Potential Number of Deaths Had There Not Been DIH Prosecutions")) +
labs(x = "Date", y = "People",
title = "Number of OD Deaths and Estimated Number of Lives Saved in Each Year",
color = "", linetype = "") +
theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank(),
panel.background = element_blank(), axis.line = element_line(colour = "black"),
legend.position = "bottom") +
scale_linetype_manual(values = c("dashed", "solid")) +
guides(color=guide_legend(nrow=2,byrow=TRUE),
linetype=guide_legend(nrow=2,byrow=TRUE))
We now run the analysis with logistic link function using a subset of the data.
sensitivity_anlys_post_tx_model_logistic_smoothed_time_subset<-gam(cbind(round(imputed_deaths), round(num_alive))~ State +
s(Time_Period_ID, bs = 'cr', by = as.factor(Region)) +
Naloxone_Pharmacy_Yes_Redefined +
lag_num_pd_w_naloxone_yes +
Naloxone_Pharmacy_No_Redefined +
lag_num_pd_w_naloxone_no +
Medical_Marijuana_Redefined +
lag_num_pd_w_med_marijuana +
Recreational_Marijuana_Redefined +
lag_num_pd_w_rec_marijuana +
GSL_Redefined +
lag_num_pd_w_gsl +
PDMP_Redefined +
lag_num_pd_w_pdmp +
Medicaid_Expansion_Redefined +
lag_num_pd_w_medicaid +
Intervention_Redefined +
lag_num_pd_w_tx,
data = data_subset,
family = "binomial")
summary(sensitivity_anlys_post_tx_model_logistic_smoothed_time_subset)
##
## Family: binomial
## Link function: logit
##
## Formula:
## cbind(round(imputed_deaths), round(num_alive)) ~ State + s(Time_Period_ID,
## bs = "cr", by = as.factor(Region)) + Naloxone_Pharmacy_Yes_Redefined +
## lag_num_pd_w_naloxone_yes + Naloxone_Pharmacy_No_Redefined +
## lag_num_pd_w_naloxone_no + Medical_Marijuana_Redefined +
## lag_num_pd_w_med_marijuana + Recreational_Marijuana_Redefined +
## lag_num_pd_w_rec_marijuana + GSL_Redefined + lag_num_pd_w_gsl +
## PDMP_Redefined + lag_num_pd_w_pdmp + Medicaid_Expansion_Redefined +
## lag_num_pd_w_medicaid + Intervention_Redefined + lag_num_pd_w_tx
##
## Parametric coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -9.9142839 0.0169704 -584.212 < 2e-16 ***
## StateAlaska 0.4656840 0.0395133 11.786 < 2e-16 ***
## StateArizona 0.3260481 0.0181969 17.918 < 2e-16 ***
## StateArkansas -0.3816820 0.0264243 -14.444 < 2e-16 ***
## StateCalifornia 0.1374486 0.0253367 5.425 5.80e-08 ***
## StateColorado 0.2421707 0.0256020 9.459 < 2e-16 ***
## StateConnecticut 0.2023745 0.0227399 8.900 < 2e-16 ***
## StateDelaware 0.1746091 0.0329617 5.297 1.18e-07 ***
## StateFlorida 0.5171831 0.0179741 28.774 < 2e-16 ***
## StateGeorgia 0.2669132 0.0206240 12.942 < 2e-16 ***
## StateHawaii -0.3190920 0.0385663 -8.274 < 2e-16 ***
## StateIdaho -0.3220520 0.0325556 -9.892 < 2e-16 ***
## StateIllinois 0.2045205 0.0202426 10.103 < 2e-16 ***
## StateIndiana -0.1287957 0.0204247 -6.306 2.87e-10 ***
## StateIowa -0.6289116 0.0277864 -22.634 < 2e-16 ***
## StateKansas -0.2274021 0.0245331 -9.269 < 2e-16 ***
## StateKentucky 0.5383433 0.0196523 27.393 < 2e-16 ***
## StateLouisiana 0.4179674 0.0188506 22.173 < 2e-16 ***
## StateMaine 0.1200546 0.0349900 3.431 0.000601 ***
## StateMaryland -1.3997811 0.0323084 -43.326 < 2e-16 ***
## StateMassachusetts -0.0196295 0.0199764 -0.983 0.325788
## StateMichigan -0.0430901 0.0201842 -2.135 0.032774 *
## StateMinnesota -0.5445597 0.0251231 -21.676 < 2e-16 ***
## StateMississippi -0.0899853 0.0230006 -3.912 9.14e-05 ***
## StateMissouri 0.2681887 0.0207599 12.919 < 2e-16 ***
## StateMontana -0.0404765 0.0382705 -1.058 0.290219
## StateNebraska -0.8858219 0.0371221 -23.862 < 2e-16 ***
## StateNevada 0.6634038 0.0263718 25.156 < 2e-16 ***
## StateNew Hampshire 0.1264438 0.0289908 4.362 1.29e-05 ***
## StateNew Jersey 0.1877752 0.0190435 9.860 < 2e-16 ***
## StateNew Mexico 0.8956733 0.0258327 34.672 < 2e-16 ***
## StateNew York -0.2818616 0.0188755 -14.933 < 2e-16 ***
## StateNorth Carolina 0.2273335 0.0168711 13.475 < 2e-16 ***
## StateNorth Dakota -1.3175979 0.0667482 -19.740 < 2e-16 ***
## StateOhio 0.5456305 0.0183878 29.674 < 2e-16 ***
## StateOklahoma 0.3790127 0.0205461 18.447 < 2e-16 ***
## StateOregon 0.0479041 0.0289839 1.653 0.098375 .
## StatePennsylvania 0.5443358 0.0198662 27.400 < 2e-16 ***
## StateRhode Island 0.0880039 0.0325676 2.702 0.006888 **
## StateSouth Carolina 0.1159933 0.0196563 5.901 3.61e-09 ***
## StateSouth Dakota -1.0671710 0.0567210 -18.814 < 2e-16 ***
## StateTennessee 0.3589872 0.0172473 20.814 < 2e-16 ***
## StateTexas 0.0417948 0.0194628 2.147 0.031760 *
## StateUtah -0.0404964 0.0260012 -1.557 0.119356
## StateVermont -0.1079730 0.0434154 -2.487 0.012884 *
## StateVirginia -0.1214341 0.0194204 -6.253 4.03e-10 ***
## StateWashington 0.3361612 0.0252409 13.318 < 2e-16 ***
## StateWest Virginia 0.6700384 0.0222475 30.117 < 2e-16 ***
## StateWisconsin 0.0452589 0.0202226 2.238 0.025219 *
## StateWyoming 0.0198349 0.0416322 0.476 0.633768
## Naloxone_Pharmacy_Yes_Redefined -0.0451257 0.0135442 -3.332 0.000863 ***
## lag_num_pd_w_naloxone_yes -0.0348725 0.0036484 -9.558 < 2e-16 ***
## Naloxone_Pharmacy_No_Redefined 0.1033996 0.0105883 9.765 < 2e-16 ***
## lag_num_pd_w_naloxone_no -0.0142420 0.0012446 -11.443 < 2e-16 ***
## Medical_Marijuana_Redefined 0.0208996 0.0101291 2.063 0.039081 *
## lag_num_pd_w_med_marijuana -0.0050364 0.0010727 -4.695 2.66e-06 ***
## Recreational_Marijuana_Redefined -0.1575028 0.0385241 -4.088 4.34e-05 ***
## lag_num_pd_w_rec_marijuana -0.0170023 0.0185896 -0.915 0.360394
## GSL_Redefined -0.0262166 0.0102647 -2.554 0.010648 *
## lag_num_pd_w_gsl 0.0005865 0.0025692 0.228 0.819417
## PDMP_Redefined -0.0457857 0.0078718 -5.816 6.01e-09 ***
## lag_num_pd_w_pdmp 0.0074035 0.0008069 9.175 < 2e-16 ***
## Medicaid_Expansion_Redefined 0.0467635 0.0104219 4.487 7.22e-06 ***
## lag_num_pd_w_medicaid 0.0173406 0.0029639 5.851 4.90e-09 ***
## Intervention_Redefined 0.0151905 0.0070978 2.140 0.032341 *
## lag_num_pd_w_tx -0.0180270 0.0006636 -27.165 < 2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Approximate significance of smooth terms:
## edf Ref.df Chi.sq p-value
## s(Time_Period_ID):as.factor(Region)Midwest 8.530 8.932 3285 <2e-16 ***
## s(Time_Period_ID):as.factor(Region)Northeast 8.641 8.960 2127 <2e-16 ***
## s(Time_Period_ID):as.factor(Region)South 8.681 8.969 3411 <2e-16 ***
## s(Time_Period_ID):as.factor(Region)West 8.911 8.998 1468 <2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## R-sq.(adj) = 0.882 Deviance explained = 87.3%
## UBRE = 5.5488 Scale est. = 1 n = 1480
#compute the full dataset including basis functions
full_df_w_basis_functions_sensitivity_anlys_post_tx_logistic_smoothed_time_subset <-
data.frame(predict(sensitivity_anlys_post_tx_model_logistic_smoothed_time_subset, type = "lpmatrix"))
#estimate the 95% CI and S#
coefficient_values_sensitivity_anlys_post_tx_logistic_smoothed_time_subset <- coef(sensitivity_anlys_post_tx_model_logistic_smoothed_time_subset)
sensitivity_anlys_post_tx_sd_and_ci_logistic_smoothed_time_subset_w_cov <-
compute_sd_and_CI_logistic_link(as.matrix(full_df_w_basis_functions_sensitivity_anlys_post_tx_logistic_smoothed_time_subset),
(data_subset$prop_dead),
coefficient_values_sensitivity_anlys_post_tx_logistic_smoothed_time_subset,
d = ncol(full_df_w_basis_functions_sensitivity_anlys_post_tx_logistic_smoothed_time_subset),
return_full_cov = TRUE)
#create a data set for plot
sensitivity_anlys_post_tx_sd_and_ci_logistic_smoothed_time_subset <-
sensitivity_anlys_post_tx_sd_and_ci_logistic_smoothed_time_subset_w_cov[[1]]
colnames(sensitivity_anlys_post_tx_sd_and_ci_logistic_smoothed_time_subset) <- c("conf.low", "estimate", "conf.high",
"exp_lb", "exp_coef", "exp_ub", "sd")
sensitivity_anlys_post_tx_sd_and_ci_logistic_smoothed_time_subset$term <- rownames(sensitivity_anlys_post_tx_sd_and_ci_logistic_smoothed_time_subset)
sensitivity_anlys_post_tx_sd_and_ci_logistic_smoothed_time_subset$ci_95 <-
paste("95% CI = (", format(round(sensitivity_anlys_post_tx_sd_and_ci_logistic_smoothed_time_subset$conf.low, 3), nsmall = 3), ", ",
format(round(sensitivity_anlys_post_tx_sd_and_ci_logistic_smoothed_time_subset$conf.high, 3), nsmall = 3), ")", sep = "")
dwplot(sensitivity_anlys_post_tx_sd_and_ci_logistic_smoothed_time_subset[51:66,], colour = "black") +
theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank(),
panel.background = element_blank(), axis.line = element_line(colour = "black")) +
geom_vline(aes(xintercept = 0), linetype = "dashed") +
labs(y = "Term", x = "Coefficients and 95% Confidence Intervals",
title = "Coefficient of Analysis With Smoothed Time Effects,
Linear Intervention, Subset Data") +
scale_color_grey() +
geom_text(sensitivity_anlys_post_tx_sd_and_ci_logistic_smoothed_time_subset[51:66,],
mapping = aes(label = format(round(estimate, 3), nsmall = 3), x = 0.55, y = 16:1), size = 3) +
geom_text(sensitivity_anlys_post_tx_sd_and_ci_logistic_smoothed_time_subset[51:66,],
mapping = aes(label = ci_95, x = 0.9, y = 16:1), size = 3) +
xlim(-.5, 1.1)
table_of_RR_lin_eff_logistic_all_yr_subset <- sensitivity_anlys_post_tx_sd_and_ci_logistic_smoothed_time_subset[51:66,]
table_of_RR_lin_eff_logistic_all_yr_subset$estimate <- round(exp(table_of_RR_lin_eff_logistic_all_yr_subset$estimate), 3)
table_of_RR_lin_eff_logistic_all_yr_subset$conf.low <- exp(table_of_RR_lin_eff_logistic_all_yr_subset$conf.low)
table_of_RR_lin_eff_logistic_all_yr_subset$conf.high <- exp(table_of_RR_lin_eff_logistic_all_yr_subset$conf.high)
table_of_RR_lin_eff_logistic_all_yr_subset$ci_95 <- paste("95% CI = (",
format(round(table_of_RR_lin_eff_logistic_all_yr_subset$conf.low, 3),
nsmall = 3), ", ",
format(round(table_of_RR_lin_eff_logistic_all_yr_subset$conf.high, 3),
nsmall = 3), ")", sep = "")
# write.csv(table_of_RR_lin_eff_logistic_all_yr_subset, "./Data/table_of_RR_lin_eff_logistic_all_yr_subset_5_13_22.csv")
table_of_RR_lin_eff_logistic_all_yr_subset
## conf.low estimate conf.high exp_lb
## Naloxone_Pharmacy_Yes_Redefined 0.7975356 0.956 1.1456561 0.7975356
## lag_num_pd_w_naloxone_yes 0.9191672 0.966 1.0146484 0.9191672
## Naloxone_Pharmacy_No_Redefined 0.8671666 1.109 1.4181075 0.8671666
## lag_num_pd_w_naloxone_no 0.9707764 0.986 1.0011758 0.9707764
## Medical_Marijuana_Redefined 0.8404603 1.021 1.2406119 0.8404603
## lag_num_pd_w_med_marijuana 0.9782338 0.995 1.0120053 0.9782338
## Recreational_Marijuana_Redefined 0.6278541 0.854 1.1623476 0.6278541
## lag_num_pd_w_rec_marijuana 0.8793360 0.983 1.0992011 0.8793360
## GSL_Redefined 0.8431123 0.974 1.1254939 0.8431123
## lag_num_pd_w_gsl 0.9761370 1.001 1.0256488 0.9761370
## PDMP_Redefined 0.8120778 0.955 1.1236560 0.8120778
## lag_num_pd_w_pdmp 0.9860868 1.007 1.0292370 0.9860868
## Medicaid_Expansion_Redefined 0.9032576 1.048 1.2156447 0.9032576
## lag_num_pd_w_medicaid 0.9743360 1.017 1.0625593 0.9743360
## Intervention_Redefined 0.8799492 1.015 1.1714849 0.8799492
## lag_num_pd_w_tx 0.9677675 0.982 0.9967149 0.9677675
## exp_coef exp_ub sd
## Naloxone_Pharmacy_Yes_Redefined 0.9558774 1.1456561 0.092399585
## lag_num_pd_w_naloxone_yes 0.9657285 1.0146484 0.025211573
## Naloxone_Pharmacy_No_Redefined 1.1089344 1.4181075 0.125471269
## lag_num_pd_w_naloxone_no 0.9858590 1.0011758 0.007865877
## Medical_Marijuana_Redefined 1.0211195 1.2406119 0.099339364
## lag_num_pd_w_med_marijuana 0.9949763 1.0120053 0.008658247
## Recreational_Marijuana_Redefined 0.8542744 1.1623476 0.157114571
## lag_num_pd_w_rec_marijuana 0.9831414 1.0992011 0.056931590
## GSL_Redefined 0.9741241 1.1254939 0.073693121
## lag_num_pd_w_gsl 1.0005867 1.0256488 0.012621863
## PDMP_Redefined 0.9552466 1.1236560 0.082843573
## lag_num_pd_w_pdmp 1.0074309 1.0292370 0.010925678
## Medicaid_Expansion_Redefined 1.0478742 1.2156447 0.075770939
## lag_num_pd_w_medicaid 1.0174919 1.0625593 0.022112113
## Intervention_Redefined 1.0153064 1.1714849 0.073000832
## lag_num_pd_w_tx 0.9821346 0.9967149 0.007518587
## term
## Naloxone_Pharmacy_Yes_Redefined Naloxone_Pharmacy_Yes_Redefined
## lag_num_pd_w_naloxone_yes lag_num_pd_w_naloxone_yes
## Naloxone_Pharmacy_No_Redefined Naloxone_Pharmacy_No_Redefined
## lag_num_pd_w_naloxone_no lag_num_pd_w_naloxone_no
## Medical_Marijuana_Redefined Medical_Marijuana_Redefined
## lag_num_pd_w_med_marijuana lag_num_pd_w_med_marijuana
## Recreational_Marijuana_Redefined Recreational_Marijuana_Redefined
## lag_num_pd_w_rec_marijuana lag_num_pd_w_rec_marijuana
## GSL_Redefined GSL_Redefined
## lag_num_pd_w_gsl lag_num_pd_w_gsl
## PDMP_Redefined PDMP_Redefined
## lag_num_pd_w_pdmp lag_num_pd_w_pdmp
## Medicaid_Expansion_Redefined Medicaid_Expansion_Redefined
## lag_num_pd_w_medicaid lag_num_pd_w_medicaid
## Intervention_Redefined Intervention_Redefined
## lag_num_pd_w_tx lag_num_pd_w_tx
## ci_95
## Naloxone_Pharmacy_Yes_Redefined 95% CI = (0.798, 1.146)
## lag_num_pd_w_naloxone_yes 95% CI = (0.919, 1.015)
## Naloxone_Pharmacy_No_Redefined 95% CI = (0.867, 1.418)
## lag_num_pd_w_naloxone_no 95% CI = (0.971, 1.001)
## Medical_Marijuana_Redefined 95% CI = (0.840, 1.241)
## lag_num_pd_w_med_marijuana 95% CI = (0.978, 1.012)
## Recreational_Marijuana_Redefined 95% CI = (0.628, 1.162)
## lag_num_pd_w_rec_marijuana 95% CI = (0.879, 1.099)
## GSL_Redefined 95% CI = (0.843, 1.125)
## lag_num_pd_w_gsl 95% CI = (0.976, 1.026)
## PDMP_Redefined 95% CI = (0.812, 1.124)
## lag_num_pd_w_pdmp 95% CI = (0.986, 1.029)
## Medicaid_Expansion_Redefined 95% CI = (0.903, 1.216)
## lag_num_pd_w_medicaid 95% CI = (0.974, 1.063)
## Intervention_Redefined 95% CI = (0.880, 1.171)
## lag_num_pd_w_tx 95% CI = (0.968, 0.997)
date_data_subset <- data_subset[, c("Time_Period_ID", "Time_Period_Start")]
date_data_subset <- date_data_subset[!duplicated(date_data_subset),]
var_cov_mat_beta_logistic_tx_subset <- as.matrix(sensitivity_anlys_post_tx_sd_and_ci_logistic_smoothed_time_subset_w_cov[[2]][
c("Intervention_Redefined", "lag_num_pd_w_tx"), c("Intervention_Redefined", "lag_num_pd_w_tx")])
attr_deaths_est_logistic_smoothed_time_lin_post_subset <- attr_death_compute(
data_subset,
sensitivity_anlys_post_tx_sd_and_ci_logistic_smoothed_time_subset,
tx_name = c("Intervention_Redefined", "lag_num_pd_w_tx"),
var_cov_mat_beta_logistic_tx_subset)
attr_deaths_est_logistic_smoothed_time_lin_post_subset <- merge(attr_deaths_est_logistic_smoothed_time_lin_post_subset, date_data_subset,
by.x = "Time_Period", by.y = "Time_Period_ID")
attr_deaths_est_logistic_smoothed_time_lin_post_subset_summary <- attr_deaths_est_logistic_smoothed_time_lin_post_subset %>%
group_by(year = year(Time_Period_Start)) %>%
summarise(total_attr_deaths = sum(attr_deaths),
total_attr_deaths_lb = sum(attr_deaths_lb),
total_attr_deaths_ub = sum(attr_deaths_ub))
sum(attr_deaths_est_logistic_smoothed_time_lin_post_subset_summary$total_attr_deaths)
## [1] 65061.87
sum(attr_deaths_est_logistic_smoothed_time_lin_post_subset_summary$total_attr_deaths_lb)
## [1] -32979.83
sum(attr_deaths_est_logistic_smoothed_time_lin_post_subset_summary$total_attr_deaths_ub)
## [1] 163103.6
sum(attr_deaths_est_logistic_smoothed_time_lin_post_subset_summary$total_attr_deaths)/sum(data_subset$imputed_deaths)
## [1] 0.1716333
sum(attr_deaths_est_logistic_smoothed_time_lin_post_subset_summary$total_attr_deaths_lb)/sum(data_subset$imputed_deaths)
## [1] -0.08700085
sum(attr_deaths_est_logistic_smoothed_time_lin_post_subset_summary$total_attr_deaths_ub)/sum(data_subset$imputed_deaths)
## [1] 0.4302675
ggplot(attr_deaths_est_logistic_smoothed_time_lin_post_subset_summary, aes(x = year)) +
# geom_point(aes(y = attr_deaths)) +
geom_line(aes(y = total_attr_deaths, linetype = "Estimate")) +
# geom_point(aes(y = attr_deaths_lb)) +
geom_line(aes(y = total_attr_deaths_lb, linetype = "95% CI")) +
# geom_point(aes(y = attr_deaths_ub)) +
geom_line(aes(y = total_attr_deaths_ub, linetype = "95% CI")) +
labs(x = "Date", y = "Lives Saved",
title = "Estimated Number of Lives Saved Per Year
Using Smoothed Time Effects,
Log Probability of Drug Overdose Death, Linear Policy Effects",
linetype = "") +
theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank(),
panel.background = element_blank(), axis.line = element_line(colour = "black")) +
scale_linetype_manual(values = c("dashed", "solid"))
ggplot() +
geom_line(attr_deaths_est_logistic_smoothed_time_lin_post_summary,
mapping = aes(x = year, y = total_attr_deaths, linetype = "Estimate",
color = "Full Dataset")) +
geom_line(attr_deaths_est_logistic_smoothed_time_lin_post_summary,
mapping = aes(x = year, y = total_attr_deaths_lb, linetype = "95% CI",
color = "Full Dataset")) +
geom_line(attr_deaths_est_logistic_smoothed_time_lin_post_summary,
mapping = aes(x = year, y = total_attr_deaths_ub, linetype = "95% CI",
color = "Full Dataset")) +
geom_line(attr_deaths_est_logistic_smoothed_time_lin_post_subset_summary,
mapping = aes(x = year, y = total_attr_deaths, linetype = "Estimate",
color = "Subset Dataset")) +
geom_line(attr_deaths_est_logistic_smoothed_time_lin_post_subset_summary,
mapping = aes(x = year, y = total_attr_deaths_lb,
linetype = "95% CI",
color = "Subset Dataset")) +
geom_line(attr_deaths_est_logistic_smoothed_time_lin_post_subset_summary,
mapping = aes(x = year, y = total_attr_deaths_ub,
linetype = "95% CI",
color = "Subset Dataset")) +
labs(x = "Date", y = "Lives Saved",
title = "Estimated Number of Lives Saved Using Smoothed Time Effects,
Log Probability of Drug Overdose Death, Linear Policy Effects",
linetype = "", color = "") +
theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank(),
panel.background = element_blank(), axis.line = element_line(colour = "black")) +
scale_linetype_manual(values = c("dashed", "solid"))